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Summary: Schlieren photographs of a small scale approximately two-dimen- 
sional air jet in air above the critical pressure have shown the existence of an 
instability displaying an anti-symmetric pattern, and the associated sound field 
having several distinguishing characteristics, has been photographed. An elementary 
theory based upon the hypothesis that the acoustic energy originates from the inter- 
action of the stream disturbances with the shock waves of the flow, and connecting 
the stream disturbances with the radiated sound, is suggested. This is shown to be 
consistent with the physical dimensions of the phenomenon and predicts a sound 
field in reasonable agreement with that observed after certain simplifying 
assumptions have been made. 


1. Introduction 


During an examination of the flow of air jets in air by the Toepler-schlieren 
optical method"?, certain disturbances were noticed in the flow from circular jets. 
These disturbances, which were in addition to the general turbulence, were observed - 
at pressure ratios above the critical and appeared along the jet boundaries adjacent 
to the resultant cellular pattern, and involved a certain amount of the external fluid 
as well as the jet stream itself. Faint sound waves were also visible. 


To bring this form of disturbance into greater prominence, similar photographs 
were taken of a jet issuing from a narrow rectangular orifice, when it was found that 
it was possible to photograph, with greater clarity, the sound waves emanating from 
the stream in the neighbourhood of these disturbances, which were of a distinct 
alternate form. It is with the relationship between the disturbances and the observed 
sound waves in the approximately two-dimensional case that this paper is primarily 
concerned. 


Notation 
c speed of sound in ambient atmosphere 


d___jet depth (small dimension of rectangular exit) 
D, directionality of the fundamental sound field 


Paper received July 1952. 
{The Aeronautical Quarterly, Volume IV, February 1953] 
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directionality of the first harmonic sound field 
relative strength of the first harmonic 

wave number, 27/A 

p’ +m’ 

assumed number of sound waves 

Mach number 

assumed number of complete cycles of stream disturbance 
m in terms of 7 

number of cells between exit and source 
defined in Section 5S(ii) 

p in terms of 7 

pressure ratio 

critical value of R for sonic exit speed 
distances defined by Fig. 9 

cell length 

times defined by Fig. 7 

relative phase of harmonic 

defined by Fig. 9 

phase distance 

disturbance spacing at the effective source 
wavelength 

time differences in phase between S, and S, respectively and § 
time differences when 7, =7, 

velocity potential 


~ 


3 


~ 


= 


2. Apparatus 


The air from a small reservoir supplied by a small compressor entered a small 
settling chamber before entering a pipe, the end of which was shaped to form a 
nozzle. Tappings at the chamber enabled pressure and temperature measurements 
to be made; the departure from ambient conditions of air temperature within the 
chamber was small. 


As no air-drying plant was available, it was inevitable that the air was used in 
a comparatively humid state, and that there was little control over that state. 


Two designs of nozzle, and the pipe carrying it, were used. The total lengths 
were the same, 19 in., and the dimensions of the rectangular exit were similar, being 
0:07 in. x 0-118 in. and 0-120 in. respectively. One nozzle, designated “A ” was 
formed by swaging the end of a half-inch internal diameter drawn copper tube over 
a length of two inches. The other nozzle “ B” was formed by tapering the end of a 
one-inch square brass pipe over a length of four inches; to obtain the requisite 
pressure ratio at the exit it was necessary to reduce the larger exit dimension to 
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0-70 in. by wall insertions, which were tapered away internally over the length of 
the original contraction. 


The connections at the chamber were adapted to carry either of these nozzles. 
A fine gauze mesh close to the pipe entry was originally used as a turbulence screen; 
but as it did not appear to influence the jet flow its use was discontinued after a 
while to reduce the overall pressure losses. 


The optical system was of the conventional Toepler-schlieren type, excepi 
perhaps for the use of an auxiliary condensing lens and a small rectangular knife- 
edge window, the latter forming an effective source. Apart from eliminating 
troublesome images of the structure immediately surrounding the spark source and 
so providing a small well defined image, the knife-edge window could also be used 
to provide a control over the illumination, a method found more convenient than 
adjustment of the actual source or of the circuit components. 


The light source was provided by a spark discharge between magnesium 
electrodes, which were enclosed in a short length of “Pyrex” capillary tube to 
constrain the spark path. The flash produced was of the order of 0:2 candle-seconds 
with an effective duration of one or two microseconds. 


Two major requirements arise when the system is being used for the detection 
of very small density gradients. Firstly, the “cut-off” at the critical knife edge must 
be large, so that little light in the undisturbed state is allowed to pass, so as to 
obtain a high sensitivity. Diffraction, producing effects considered undesirable in 
this work, results in a limit to the sensitivity of the system; generally it was found 
that satisfactory results were obtained when a strip of the image about 0-007 in. wide 
was unstopped. The other requirement is that the resultant emulsion density at the 
photographic plate shall be the optimum; this was obtained by a trial and error 
adjustment of the auxiliary knife edges stopping the length of the image. Standard 
photographic technique proved satisfactory, using ortho-chromatic film and 
ordinary developing methods. 


3. Observed Jet Phenomena 


The phenomena which have been observed in the approximately two-dimen- 
sional case, where the pressure ratio exceeds the critical, may be divided into three 
groups : — 


(i) The stationary repetitive or “cellular” pattern* arising as a consequence of 
the excess pressure at the nozzle exit. 


(ii) The turbulence of the stream, including a definite alternate formation in 
addition to the general turbulence. 


*The term “cellular” has been used in order to avoid possible confusion with other regular 
or periodic characteristics, 
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(iii) The emission of sound waves of a form which enable them to be 
photographed. 


These points refer specifically to the two-dimensional case, although examina- 
tion of photographs of jets issuing from circular orifices suggests that the same type 
of phenomena. exist with some change of character, particularly as regards the 
last group. 


The nature of the observed phenomena can be seen in the photograph of Fig. 1 
(see page 109) and will now be considered in greater detail. 


' 


(i) The cellular pattern 

This characteristic is peculiar to jets issuing from an orifice at a sonic or 
supersonic velocity with an exit pressure different from that of the ambient atmos- 
phere. Considering the case of exit pressure greater than that of the ambient 
atmosphere, the flow at first expands, but as a result of the reflection at the constant 
pressure boundaries of the Mach lines, originating from the edge of the exit, the 
flow then contracts and approaches conditions very similar to those existing at the 
orifice. If it were not for turbulence setting in from the edges of the jet, the process 
would repeat indefinitely; in practice the pattern usually extends over about four 
or five cycles before becoming lost in the general turbulence. (In the three-dimen- 
sional case as many as twelve cycles have been observed.) 


25 
20 
15 
s 
d 
1:0 
1-5 
Fig. 2. 


Variation of cell spacing with excess pressure parameter. 


Owing to the depth of field observed, the density gradients associated with this 
cellular pattern are too great to be depicted by a continuous variation of illumination, 
as the optical system had been adjusted for the recording of very small gradients. 
Consequently these gradients result in either a complete loss of light or a local over- 
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exposure of the photographic plate. For this reason, any shock waves which may 
develop.in the flow are not shown; they would occur at the end of the cells growing 
inwards in the upstream direction with increasing pressure ratio, being strongest at 
the boundary”. 


The spacing of the cellular pattern, or “cell” length (Fig. 2), follows an empiri- 
cal law of a nature similar to that for jets issuing from a circular orifice. The 
relation for the cell adjacent to the exit is 


s_celllength |. _RY 


where “jet depth” means the small dimension of the rectangular exit. R denotes 
the pressure ratio, R. being the critical value to produce a sonic exit speed. (It should 
be noted that the values of pressure ratio are approximate and are as derived from 
pitot measurements at the exit, using hypodermic tubing). 


In some cases condensation shocks have appeared close to the exit; the result 
on the flow appears to be a recovery of pressure large enough to displace the normal 
cellular pattern downstream, measurements of the lengths of individual cells indi- 
cating that they have not undergone any appreciable modification. The point 
corresponding to the lowest pressure ratio is such a case. 


(ii) The turbulence of the jet 


The jet boundary close to the exit appears to be similar to those known to be 
laminar, although the exit Reynolds number of about 10° and the length of the 
approach tube, with unsteady entry conditions, indicate that the flow as a whole - 
would be turbulent. 


At the end of this laminar region, turbulence apparently of a random nature 
begins to develop steadily with increasing distance downstream to spread over a 
total angle upwards of 20°. : It is difficult to.assign a mean line to the boundaries, 
as the distance observed is not very large in comparison with the magnitude of the 
irregularities, so that the measured angles are not representative of the flow as a 
whole. 


An outstanding feature of the turbulence is the pattern of strong density 
gradients visible in the region where the turbulence proper sets in. This pattern 
involves both the stream itself, to a considerable depth and sometimes appearing to 
displace it laterally, and also a certain amount of the fluid external to the apparent 
jet boundary. The pattern takes up a very definite alternate formation and displays 
three or four periods, apparently being damped out in the developing turbulent flow. 
The introduction of “eddy viscosity” may be responsible for this, although it is 
interesting to note that Brown has observed sensitive jet disturbances which are 
damped after a period of amplification, the flow being laminar throughout. In the 


107 


OEE 
be 
a- 
pe 
he 

or 

nt 
nt 

he 
e 

SS 
r 


ALAN POWELL 


profile is such that damping takes place for that particular wavelength. The spacing 
of the pattern, both axially and laterally, increases somewhat in the downstream 
direction. 


7:0 


V 


latter case the disturbances have moved to an area in the flow where the velocity | 
| 
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Fig. 4. 
Variation of wavelength of upstream radiation. 


(iii) 


The observed sound field 

The density gradients of some of the sound waves emanating from the jet stream 
are sufficiently large to be photographed by the schlieren method under certain 
conditions. 


Most prominent of the sound waves are those passing in the upstream direction 
on either side of the jet axis, and these are in anti-phase on opposite sides of the axis 
(Fig. 1). Reference to a case where the knife edges were inclined at 45° to the jet 
axis (Fig. 6—p. 110) reveals that the intensity of these sound waves begins to fall 
away quite rapidly for angles greater than about 45° to the upstream jet axis, and 
little is detectable at angles greater than 70° or so. In some cases there is an 
indication of some radiation at this frequency passing in the downstream direction, 
but this is always weak in comparison with the upstream radiation. 


The wavelength of the upstream radiation is very dependent upon the pressure 
ratio, and is in fact closely linked with the square root of the super-pressure, as is 
the cell length. A set of experimental points is plotted in Fig. 4 and from this it 
follows that the wavelength is given to a fair approximation by the empirical relation | 
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Fig. 1. A typical photograph (R = 2°88). 


Fig. 3. Addition of reflector plates (R =2°56). 
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Fig. 5. Addition of reflector plates (R =3-17). 


Fig. 6. Photograph with knife edges inclined at 45° (top left to bottom right). (R=3-17). 
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The lower point, it is to be noted, was obtained from a photograph displaying a 
strong condensation shock. 


To observe any sufficiently intense radiation in directions normal to the jet 
axis, reflector plates were introduced so that the density gradients of such radiation 
would be detectable using a vertical knife-edge system as in Fig. 5. This arrange- 
ment showed very clearly that there was a powerful radiation originally normal to 
the jet stream. In some cases this radiation took a complicated form, as in Fig. 5, 
but in some others (for example Fig. 3) the radiation had a very simple form. It is 
worth noting that the wavelength of the normal radiation, where it is possible to 
make a fairly accurate estimate of it, is precisely half the upstream wavelength and, 
in the other rather more complex cases, it can be seen that the pattern could be 
composed by a series of such waves. 


There is also some additional radiation present in the upstream direction. It is 
of a much higher frequency, if indeed a frequency as such can be assigned to it. 
It appears to emanate from approximately the same region of the stream as the 
other radiation. 


If certain details of the geometry of the optical system are known the intensity 
of the radiated sound may be deduced from the schlieren photographs*, provided 
that a reasonable approximation can be made regarding the spatial distribution of 
the source. In the cases under consideration the assumption of a point source is a. 
sufficiently good approximation for points several wavelengths distant from the 
source. By considering the density gradients necessary to produce the observed 
contrast at the photographic plate, the intensity of the sound can be deduced. In 
this way, the sound level of the waves passing in the upstream direction in Fig. 1 is 
found as approximately 165 decibels—a powerful emission. 


4. A Suggested Mechanism 


The production of noise from the jets described in this paper has a significant 
difference in comparison with normal aerodynamic noise which has a continuous 
spectrum; this is the existence of dominant frequencies for given jet conditions. In 
view of the fact that the radiated wavelengths change from a minimum to a maxi- 
mum for a comparatively small increase in Reynolds number, and that the radiated 
wavelengths are roughly proportional to the spacing of the cellular pattern, it would 
seem that the mechanism responsible for the sound is stabilised in some way at a 
certain frequency. The following mechanism is suggested as fulfilling this require- 
ment, and at the same time providing an explanation of the characteristic acoustic 


*A note concerning this method is being prepared. 
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field. It is assumed that the sound waves moving in the upstream direction adjacent 
to the jet are of sufficient strength to govern the stability of the boundary of the 
stream close to the exit, so giving rise to embryo disturbances. These periodic 
stream disturbances, which will be carried downstream at a certain velocity, are 
assumed to be amplified in passing downstream and, after a certain distance, 
produce significant acoustic energy on traversing the shock waves associated with 
the cellular system; part of this energy will pass in the upstream direction so 
perpetuating the periodic process. 


The association of the production of acoustic energy with an instability of shock 
waves has been known for some time. Hilton”) observed that the shock waves 
on an aerofoil oscillate with the same frequency as that of the sound attributed 
to them. A similar state of affairs may have existed in some more recent experi- 
ments by Holder and North. 


It is not difficult to infer from a physical point of view in an elementary way 
that a disturbance traversing a shock wave is likely to produce acoustic energy, for 
changes in velocity or pressure on the upstream side of the shock produce variations 
in both velocity and pressure on the downstream side (as well as a movement of 
the shock position) and these may be in part relieved by radiation to a distance. 


The influence of sound on jet flow has been the subject of much experiment, 
although much remains to be done in this field. It has been shown that the nature 
of the flow can be greatly modified if sound waves are permitted to impinge upon 
the jet in the immediate neighbourhood of the exit. This is the basis of the pheno- 
mena known as sensitive jets’ or flames’, when a very distinct vortex pattern may 
be produced in the flow. These phenomena occur, of course, only under certain 
conditions, the principal one being that the flow from the exit must be laminar, 
although the critical Reynolds number may be exceeded somewhat if stream 
disturbances are carefully guarded against’. However, it has been observed that 
some influence on an initially turbulent stream is produced if the sound is sufficiently 
intense. The associated phenomenon of “edge tones” has also been observed 
using jets in a very turbulent state”). In the related case of flow past a fixed 
boundary, it has been shown that sound waves may give rise to unstable oscillations 
in a laminar boundary layer, enhancing the onset of turbulence®. The stability of 
disturbances, alternately disposed on the boundaries of a jet was first discussed 
by Savic’ and the work has been extended more recently to the compressible case 
by Pai"°’, but these works are limited to the case where the flow is wholly laminar 
and initially steady. 


Thus there is some foundation underlying the assumed energy linkages from 
stream disturbances to acoustic waves and vice versa. The next step is to show 
that the physical dimensions of the system are consistent with this mechanism and 
that it is capable of explaining the predominant features of the sound field, i.e. its 
directionality. 


112 


NOISE FROM AIR JETS 


5. On the Dimensions of the System 


If the mechanism of sound production is as assumed, the net output will be the 
sum of the outputs from more than one point, for an examination of the photo- 
graphs will show that the periodic stream disturbance system and cellular pattern 
have an appreciable strength and are coincident for some length of the stream. 
Clearly, since these individual sources are correlated in some way, they can be 
combined to form an effective source acting at a certain point. 


Consider the mechanism in greater detail. The basic conception is that a 
disturbance effectively creates a sound wave on traversing that section which is the 
location of the effective source; that this sound waves proceeds upstream and creates 
a stream disturbance in the immediate neighbourhood of the jet exit; and that this 
disturbance in turn produces a sound wave on reaching the effective source. The 
process is thus self-maintained. But several questions arise, showing that the 
mechanism as just outlined is over-simplified :— 


(i) Since both the stream disturbance and sound wave are continuous, what 
exactly is meant when it is said that a disturbance gives rise to a sound wave? 
Clearly it is possible to say that a certain phase of disturbance will be associated 
with a certain phase of sound wave, and this is sufficient for the present argument, 
in which “disturbance” and “wave” are taken to refer to a particular pair of 
associated phases. 


(ii) The sound wave, of particular phase, is to give rise to a stream disturbance 
at the origin; but this stream disturbance is not necessarily of the same phase as 
that which was responsible for the production of that particular phase of sound 
wave. This fact can be allowed for in the ensuing analysis by assuming that the 
embryo stream disturbance is lagging a certain amount “ p” behind the associated 
phase as defined. 


(iii) There is not necessarily only a single cycle of the stream disturbance 
between the end points of the mechanism. The physical dimensions and velocities 
may be such that there may be several cycles present, so that several sound waves 
are produced during the time taken for the disturbance to reach the effective source. 
Thus it will be assumed that there are “m” such complete cycles of stream dis- 
turbance between the end points. 


(iv) The same is true of the sound waves; hence “/” sound waves are assumed 
to be involved. 


(v) Acomplication arises as to whether the velocity of the disturbance remains 
constant over the distance concerned. Clearly this is not so, since the spacing of the 
alternate disturbance pattern has been observed to increase in the downstream 
direction, indicating an increasing velocity in that direction. The disturbance 
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velocity, before the disturbances are detectable, may be appreciably different from 
the later values. This variable velocity can be catered for by reference to the dis- 
turbance spacing evident in the region of the effective source. Let this spacing be 1. 


Thus, for example, the m cycles of (iii) will occupy a distance m’n where m’ is 
less than m and in general will not be an integer. 


(vi) A further definition will prove useful; the effective source is located some 
distance downstream and is associated with the spacing of the cellular pattern. Let 
it be situated a distance n’s from the orifice, where n’ is not necessarily an integer. 


© M 
N Ss 
POSITION AT TIME t, 
Wi 
AN | 
Y, 


N Ss 
POSITION AT TIME t, 
Fig. 7. 


Let the disposition of the elements of the process be as in Fig. 7 at a certain 
time ¢,, when the associated phases of sound wave (W,) and disturbance (E,) are 
about to leave S, the effective source. The stream disturbance moves downstream 
at a velocity Mc where c is the velocity of sound in the ambient atmosphere. 


Now consider the system at a slightly later time t,, when the next sound wave 
arrives at the nozzle exit N after W, left S. By the foregoing reasoning this wave 
will be a distance JA upstream from W, if A is the wavelength. Also this wave 
produce a stream disturbance of a certain phase, distant p’n from the nearest point 
in the stream having the same phase as E,, and distant m’n from it. 


114 


NOISE FROM AIR JETS 


During the time (¢,—1t,), E, will have moved to a close approximation 


n’s—IXr 
é 


Mc (t, —t,)=Mc 


so that the distance NS can be written 
n’s=(p’ +m’) — M (n’s—IA) 
or n’s=k’n—M (n’s—Id) ‘ (3) 
where k’=p’+m’. 


By this approach the three unknown variables p, m and the variation of M 
can be combined in the single parameter k’. 


In accordance with the basic assumption that sound is created as the stream 
disturbances traverse the cellular system, and therefore the effective source, a con- 
sideration of events then leads directly to the relation 


or simply that 
n=Mdr ‘ i ‘ (4) 


Thus if the assumed system is to exist it must fulfil the equations (3) and (4). 
The present state of knowledge does not make it possible to predict k’, or any of its 
components to any degree of accuracy. However, k’ can be found experimentally 
for any given case and if the assumption is made that the variations of the speed of 
the disturbance system over the distance NS can be considered similar for various 
jet conditions, then that value of k’ will be fixed unless m’ changes. By the nature 
of the system, if m’ does change, it will do so by approximately an integer indicating 
the introduction, or removal, of a complete cycle of the stream disturbance between 
N and S. 


In choosing a photograph of the “A”: series considered to be reliable and 
suitable for such measurements the following values were obtained by direct 
measurement 

n’s=0°619 in., »=0°331 in., A=0-617 in. 


Hence /=1 since 4<n’s and these, using the previous equations, result in 

In the derivation of these values, the effective source was taken as at the end 
of the third cell from the exit, this coinciding with the centres of the arcs formed 
by the sound waves. Similarly for nozzle B, where the effective source appears to 
be about 44 cells from the exit, 

n’s=1:05 in.,  1=0°465in., A=0-650 in. 
giving i=1, k’=2:87 and M =0:650. 
115 


Mc’ c 


-ALAN POWELL 


It is interesting to note that the value of k’ for the two series differs by almost 
an integer; this indicates the presence of an additional cycle of the stream distur- 
bance being carried in flow between the end points of the mechanism, since it must 
be m that has increased and not p. An examination of the photographs shows this 
to be reasonable. 


In general the eddy spacing » is the most difficult measurable quantity to esti- 
mate, so that it is perhaps desirable to compare estimated and measured values of 
this rather than A, say. From (3) and (4) 


An's 


the values of which are given in Table I for various conditions. 


Some difficulty was encountered in performing some of the measurements; in 
general all have been repeated at least twice, being repeated again if necessary, and 
a mean chosen. If account is taken of the probable accuracy of measurement (say 
+5 per cent) the agreement as shown in Table I and in Fig. 8 (for the “ A ” series) 
is as good as could be expected, especially as variations of k’ within each series 
have been ignored. 


It is interesting to note that in the case /(R—R,.)=0°592 of series “A” a 
strong condensation shock was present. The upper set of figures shows a consider- 
ably larger error when the displacement of the shock system downstream is not 
allowed for than when it is taken into account, whereas the cellular pattern itself 
displays little change (see also Fig. 8). This supports the assumption that the dis- 
turbances originate in the immediate neighbourhood of the exit, as with sensitive jets. 


TABLE I 
COMPARISON OF ESTIMATED AND MEASURED DISTURBANCE SPACING 
n’s Error: 
d d M d per cent. 
6-26 6:35 0-528 3-45 3:36 
1-133 5:97 5:97 0-546 3-26 3-18 2:6 
Nozzle 0-985 5-23 5-24 0-535 2:81 
A 0-819 4-20 4-47 0-522 2:24 2:30 +3-0 
0-592 3:36 ( 3:36] 206 {1-79} -13-2 
| 4-03 | 
1-304 6:89 10-91 0-677 4-66 4-76 4$2°1 
Nozzle 1-220 6-02 9-86 0-732 4-40 4-40 0 
B 1-070 5-41 8-75 0-651 3-87 a one 


0-900 4-60 7:25 0-650 4°31 3:17 
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Fig. 8. Variation of stream disturbance wavelength. 


6. On the Directionality of the Sound Field 

The directionality of the sound field is a feature dependent upon the nature 
of each individual source and upon the correlation between neighbouring sources. 
Each boundary of the jet is considered separately, it being assumed that each source 
acts as a “simple” source of classical acoustical theory, since there is no evidence 
to indicate that each individual source is fundamentally of a more complex nature. 
However, since the process is known only to be periodic at a certain fundamental 
frequency, it is not impertinent to assume the presence of harmonics, only the first 
of which will be considered in this analysis, since there is no experimental evidence 
to suggest the presence of higher harmonics of intensity comparable with that of the 
fundamental or first harmonic. 


The presence of three or four neighbouring acoustic sources possibly represents 
the actual conditions better than two or five such sources, as it is over such a distance 
that both the disturbance system and cellular pattern are most prominent and 
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DOWNSTREAM 


So 
Fig. 9. 


coincident, especially in those cases where the field as a whole was investigated, as 
in Fig. 5. It should be noted, however, that the assumption of a number of sources 
of equal strength is at best an approximation to the actual state of affairs. 


In order to simplify the analysis each source is considered as a point source; 
but it can be seen that this is really the same thing, so far as directionality is con- 
cerned, as a transverse line source, since this can be represented by a line-distribution 
of point sources". If the wavelength of the fundamental is 


A=2t/k 


then each source of strength 4x and 4cH for the fundamental and first harmonic 
respectively can be represented by the velocity potential 


r r 


where the velocity of sound in the ambient medium is c and r is the distance from 
the source considered. It is to be understood that the real part of the complex 
function is to be taken. 


Let the three sources be S,, S, and S (Fig. 9). As the disturbance system passes 
in the direction S, to S,, source S, will be a certain time 7, leading in phase, and 
S, a time 7, lagging in phase, relative to the source S. 


Then the velocity potential at a point P due to the three sources can be 
written as 


r r 
1 


Here z is the lag of the harmonic behind the fundamental. 
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As in the analysis the harmonic may be kept separate from the fundamental, 
due to their similarity in equation (8), it is possible to consider only the latter in 
detail and simply to deduce the corresponding expression when required. 


Referring to the figure, if r can be considered moderately large in comparison 
to s, the following approximations can be made 


r,=r+scosp 


| 
r.=r—scosB 


A simplification will be introduced in that it will be assumed that 


71) Say, 
with the abbreviation 


Then the potential for the fundamental becomes 


r{r—scos B}e™ + {r? — s* cos? 8} + r{r+scos B}e-™ 
— r) 
r (r? — s? cos? 8) 


= {a +2.cos kg) cos k (ct ~r) cos B sin kq sin k (ct -n} 


on taking the real part with (s/r)? small, so 


where [D;}?=(1+2 cos /9+ (= cos sin kq) (12) 
_, 2scos sin kq 
© — 1 
and ké,=tan r(1+2coskq)’ 


The corresponding expressions for the harmonic are 


cos 2k (ct — 2), 


with D,,° and 6, obtained from equations (12) with 2k replacing k in those equations. 


The symbols D,’ and D,,* define the directionality of the fundamentals and first 
harmonic due to the three related sources acting together. The distances 6 are the 
amounts by which the combined potential is caused to lag. Hence as s/r is small, 
we have approximately 


Dé= = + 3 k (cr, —s cos f) 
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and ; D,?= = + cos 2k (cr, — s cos 


— 


which in virtue of (3) become 


9 
2 cos (1 cos 8) 
(13) 
D,? = + £ cos —*(1 —M cos f) 
> " 


In a similar manner the directional properties of four correlated sources are 
found to be 


D,*=cos (1—M cos £)cos (1—Mcos 
(14) 


D,*= cos—~(1 cos /3) cos —Mcos 


If typical values are introduced into these equations, the theoretical direction- 
ality can be compared with experimental observations. Fig. 10 shows the estimated 
field corresponding to the case of Fig. 5, for which s/7(observed)=0-538 and 
M=0-604. Since the eddy systems are in anti-phase on each side of the stream, 
it follows that the radiation at fundamental frequency must follow suit. The 
directions at which the bulk of the fundamental and harmonic radiation occurs is 
in substantial agreement with experiment, the major divergence from observation 
being that the radiation from three sources directed downstream is over-estimated. 
This might be taken to suggest that the assumption of four sources is somewhat 


3 4 

J t 
Bro 
if 

B=0 


Fig. 10. 


Theoretical directionality of the fundamental D, and harmonic D,, for three and 
four correlated sourcés. 
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nearer to fact than three. The angle at which the bulk of the harmonic is directed 
is a few degrees from the normal to the jet axis; the error is comparatively small 
and is probably attributable to slight inaccuracies of the initial data. 


The foregoing analysis does not include a consideration of the downstream 
radiation of much shorter wavelength. Its origin is at present a matter of conjecture. 
The elementary theory suggested here is unlikely to produce results reliable enough 
to test its validity, owing to the susceptibility of the higher harmonics to small 
variations of the parameters involved; the latter are only known approximately. 
Alternatively, this radiation may arise from a similar mechanism involving the 
fine-grained turbulence in its early stages, or possibly of the nature discussed 
theoretically in a recent paper by Lighthill"*), whereby acoustic energy is produced 
by the interaction of turbulent velocity fluctuations with a transverse velocity 
gradient. 


7. Conclusions 


Strong sound waves have been observed to emanate from a two-dimensional 
jet at pressures above the critical. The production of sound is associated with a 
stream disturbance, displaying an alternate formation, by virtue of interaction with 
the shock waves of the stream. The hypothesis that the stream disturbances 
originate as a result of sound waves impinging upon the stream in the immediate 
region of the exit is supported by the agreement of the estimated dimensions of the 
system with the experimental results. The distribution of the sound field can be 
predicted in an approximate manner and this is in substantial agreement with the 
photographed field. 


One important result is the fact that in certain circumstances sound waves are 
capable of exerting a marked influence on an unstable flow at Reynolds numbers 
far in excess of those associated with “sensitive” jets, in which case the sound 
originates from an independent source. 


Much further work remains to be done; the placing of this qualitative theory 
on a quantitative basis, involving the estimation of noise intensity arising from 
disturbance /shock wave interaction, the precise mechanism by which sound induces 
disturbances in an unstable boundary layer and the rate of growth of such 
disturbances when once created. Sensitive jet phenomena may be helpful, particu- 
larly from the point of view of the initiation of the stream disturbances*. 


The relation to larger scale, axially symmetric jet flow should be considered, 
especially with regard to the loud “ hooting ” observed") to emanate from such jets 
above the critical pressure. Experimental work investigating this has been 
initiated*. 


*Note added in proof: See the author’s papers ““On Edge Tones and Associated Phenomena ” 
and “ The Mechanism and Reduction of Choked Jet Noise ” 
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Some Aspects of Laminar Boundary Layer 
Separation in Compressible Flow 
with no Heat Transfer 
to the Wall 


G. E. GADD, B.A. 
(Beit Scientific Research Fellow) 


Summary: An analysis has been made which suggests that, with the types of 
pressure distribution most usual in practice and free stream Mach numbers up to 10, 
no serious errors would be introduced into the calculation of the laminar separation 
point by the assumption that o, the Prandtl number, and w, the index of variation 
of viscosity with absolute temperature, are equal to unity. (Typical actual values of 
o and » for air are 0°72 and 0-89 respectively). 


1. Introduction 

The laminar compressible. boundary layer equations become considerably 
simplified when two parameters, i.e. 7, the Prandtl number »C,/k, and w, the index 
of variation of viscosity with absolute temperature (u Oc 0”), are assumed equal to 
unity. Having made these simplifying assumptions Stewartson and Illingworth”: * 
obtained independently a transformation between a compressible boundary layer 
with zero heat transfer to the wall and with any pressure distribution and an 
incompressible boundary layer with a related pressure distribution. The 
transformation enables conditions in the compressible case to be estimated, as any 
incompressible case can be solved approximately. 


Since, however, typical values of « and w for air are 0-72 and 0-89 respectively, 
the Stewartson transformation cannot apply accurately to the practical case. The 
aim of the following theoretical investigation is to determine the effects which fairly 
small departures of o and » from unity have on the separation distances of boundary 
layers with no heat transfer and given adverse pressure distributions. 


(Such work is likely to be useful more for the qualitative understanding than 
the quantitative prediction of boundary layer flows. In practice it will usually be 
impossible to predict adverse pressure distributions, as these will usually only occur 
because of shock wave-boundary layer inter-actions). 


Firstly the Stewartson transformation and another exact solution, that for the 
artificial case o=0, are considered. The possibility of obtaining further exact 
solutions on the same lines is discussed and judged to be remote. 


Paper received September 1951. 
[The Aeronautical Quarterly, Volume IV, February 1953] 
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Next an approximate method is developed for the case with w=1 and o near 
unity, and finally the case with o=1 and © near unity is. dealt with 
by a similar method. 


The approximate methods developed are applicable with any external velocity 
distribution for which the leading edge is not a stagnation point, but in all 
computations-a constant adverse gradient is assumed. (With this distribution the 
accurate solution obtained by Howarth for the incompressible constant adverse 
gradient case forms a convenient standard of comparison). Results indicate that 
even for free stream Mach numbers as large as 10 the change in separation distance 
as compared with the e=w=1 case is probably less than 40 per cent. when o=1 
and 0:7 <o <1, and less than 25 per cent. when o=1 and 0-7< 0 <1. 


These conclusions receive support from Young’s approximate method" based 
on the momentum equation. 


They seem to indicate that with the types of pressure distribution which usually 
cause separation in practice (where the gradient suddenly becomes adverse and 
separation occurs shortly afterwards), the theoretical position of the laminar 
separation point, with free stream Mach numbers up to 10 and zero heat transfer, 
would be approximately the same whatever values within the range 0-7-1-0 were 
assumed for o and w. 


Notation 


Suffixes: Suffix “d” refers to datum (main stream) conditions (normally at the 
leading edge of the boundary layer). Suffix 1 refers to the external flow (where 
n=1), except for 7, (which should not cause confusion, since > is zero when n=1), 
and x... Suffix “w” refers to wall conditions. 


k coefficient of heat conduction 

p pressure 
t a parameter proportional to viscous stress (see equations (17) and 

just above them) 

u___- velocity component in x-direction 

v velocity component in y-direction 

x distance measured along surface 

xX, separation distance 
Xa separation distance with w=oc=1 

y distance measured normal to surface 
C, specific heat at constant pressure 
C, specific heat at constant volume 

enthalpy (=JC,9) 


J mechanical equivalent of heat 
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M number 
T viscous stress at x in the »=o=1, u,(x) boundary layer 


ee defined by equations (14) 


x 
Y C,/C, 
n 


© absolute temperature 

coefficient of viscosity 

pal’ (with w= 1) 

p/p 

€ parameter defining shape of boundary layer profile in Howarth 
approximate method (see equations (17) and just above them) 

o Prandtl number pC,/k : constant 

p density 

7 viscous stress 

7 fora c=w=1, u(x) boundary layer 

x distance measured along surface for a *=«=1 boundary layer such 
that the external velocity at y=u,(x) 

defined by » 0c 9": assumed constant 


Terms not defined in the Notation will be found defined close to where they 
appear in the paper. 


1.1. The governing equations 


The equations of motion, continuity, and energy for the steady boundary layer 
in compressible flow are 


ou Ou _ 0 
pu + pv 


Ox dy dx (1) 
dy 
0: 


where / is the enthalpy, JC,9. 
With the independent variables changed from x, y to x, u these reduce to 
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al ar al 
and (l- 5, t +? +u (5) 
where T= p ay), 
the viscous stress. This transformation is due to Crocco™). 
The boundary conditions for equation (4) are 
u=u,7=0 ; u=0, 


where suffix 1 refers to the external flow and suffix “w ” to the wall. (The condition 
at the wall follows from equation (1) which shows that 07/d0y=dp/dx there). 


The boundary condition for equation (5) at u=u, is 1=1,. The wall boundary 
condition depends on the conditions imposed; for no heat transfer to the wall it is 


and the relation between enthalpy and local Mach number is 
u? 
(y-- 1) M? . . (7) 


It is also assumed that viscosity varies with enthalpy according to the relation 


= (7) © constant (8) 
In the external flow 
2 2 
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Finally two useful relations which follow from the foregoing are 


dx Pall, dx 1+ — pall, ax (11) 


2. Exact Solutions 
2.1 Case with c=o=1. 

By a transformation due to Stewartson and Illingworth the case with s=o=1 
can be reduced to that of an incompressible boundary layer with a related external 
velocity distribution, as follows. 


When o=1 and there is no heat transfer to the wall the solution of the enthalpy 
equation (5) is 


1+ 


so it follows from equations (7) to (12) that the equation of motion plus continuity 
(4) becomes in terms of variables x, n=u/u, 


( dp ( y-1 Or 
T [1 +( M, + - pat, (1 — 1+ —— 


1, ar od 
— Fil +[o (y—1)+y] Ma? pa dy 4 


with the boundary conditions 


7=0 n=0, P ita ty (14% 


When w=1 as well as o this equation can be transformed by the Stewartson 
transformation substitutions 


u,=U, hy | 


3y-1 
dx = Gy 
dx 
2y—1 
) 
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077 au, OF; 3 OF 
to — Pa Pa ay n*) an "Pa ay =0 
with the boundary conditions 
07; dU , 
; n=0, PapaU ,” 


Hence 7; is the viscous stress at (X, U) in an incompressible boundary layer, 
where U is the X component velocity at U/U,=n, where U,(X) is the external 
velocity distribution, and where the datum density and viscosity are the same as 
for the compressible boundary layer. 


This related incompressible case can be solved by one of the approximate 
methods available, such as Howarth’s* or Thwaites’. 


Hence the position of separation for any compressible boundary layer with 
o==1 and no heat transfer can be determined, since separation occurs when 7, 
and therefore 7;, is zero at the wall. It is found that with an external velocity 
distribution of the form u,=ua— 8x, 8 constant, the value of u,/ua at separation 
increases from 0-880 when M,=0 to quite near 1:0 at large Mach numbers (see 
Table I, which was computed by the method of Thwaites). 


2.2. Case with «=0 

This is an entirely artificial case because not only is there no real fluid with 
o«=0, but even if there were, the boundary layer equations as used here would not 
be applicable, since the temperature boundary layer would not be of the same order 
of thickness as the velocity boundary layer. Nevertheless the following solution is 
of some interest because it exemplifies an exact method and also, it helps to indicate 
the trend of variation of solutions of the boundary layer equations with o. 


When o=0 the enthalpy equation (5) becomes 


aa , 9. 2 (2,) 
Oudu our 


so that 0J/du=O0 for no heat transfer, i.e. 


1=1,=1| 1 + (1 me |. 
d 


It follows from equations (11) and (12) that the substitutions 
‘Jig, A= 24 
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reduce the equation of motion plus continuity (4) to the form 


0*g du, 0g 0g 


=0 


with the boundary conditions, 


0g du 

This is the same equation as would result from the substitution 


7 (X’, u) 
V (papa) 


2(X’,u)= 


where 7; is the viscous stress in an incompressible boundary layer with external 
velocity distribution u, (X’). Thus the problem has been reduced to one which can 
be solved by approximate methods. 


The solution shows that, with high Mach numbers, separation (which occurs 
when 7, and hence g, is zero at the wall), takes place considerably downstream of its 
position with low Mach numbers, because the pressure gradients in the related 
incompressible boundary layer are made progressively less adverse with increase of 
Mach number. Results calculated by the method of Thwaites for the case with 
u, =Ua— Bx, B constant, and »=1, appear in Table I. Comparison of this solution 
with that for c=w=1 suggests that the trend of variation of separation distance with . 
o is for the distance to increase with decreasing 7, a conclusion which agrees with 
that of Section 3. 


2.3. Possibility of further exact solutions 
With o=1 and no heat transfer to the wall it might be thought that, since in 
the equation of motion plus continuity (13) the only unknown is 7, a more general 


TABLE I 


“1 aT SEPARATION AGAINST M, FOR CONSTANT EXTERNAL VELOCITY 
u 


d 
GRADIENT, No HEAT TRANSFER CASES WITH © =o7=1 AND wo=1, o=0. 


M, 0 (Small) 2 4 6 8 10 


o=1 | 0.880+0.0092M. 0.910 0.940 0.957 0.969 0.977 
| 0.880—0.0066 0.854 0.782 0.734 0.712 0.704 
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form of the Stewartson transformation could perhaps be developed which would 
reduce this equation to the incompressible form for values of » other than unity. 
This is probably not practicable, as the following argument shows. 


The most general transformation would consist in writing f(x,»)7, for 7 and 
F (x)U, for u, and changing the variables »,x to N,X. This makes 


(2%) 


so that unless X, (0N/0n),, and f are functions of x only the term of equation (13) 
involving 7? (0?z/07*) will transform to several terms of which all except one will 
differ in form from any of the terms in the incompressible equation. Moreover, 
they will be unable to cancel with any contributions from the other terms of equation 
(13). Hence it does not seem possible to reduce this equation to the incompressible 
form unless X, (ON/0»),, and f are independent of 7. But if (0N/d»), is 
independent of 7, N must be equal to 7 since obviously N must be 0 when »=0, and 
also, to satisfy the outer edge boundary condition for the transformed equation, 


N must be 1 when »=1. Hence the method appears to be restricted to changing 
the x variable to X, where X is a function of x only, and writing u,=—F(x)U,, 
(x) 7. 


With these restrictions it would seem to be impossible to transform away the 
factor [1 + (1 — ?u,?/ua")4 1) M,?]'~° in equation (13) if » #1, so that in general 
the equation cannot apparently be reduced to the incompressible form. It is to be 
expected, therefore, that cases with ~=1 but » 1 will be susceptible only to 
approximate, or to possibly arduous numerical methods of solution. 


If « ~1 the situation in general appears to be even worse, for then, except for 
the artificial case «=0 already dealt with, there is no obvious solution to the 
enthalpy equation (5). Hence here too it is doubtful whether convenient exact 
solutions can be found. 


Accordingly in what follows approximate methods of solution are considered. 


3. Approximate Solution for Boundary Layers with no Heat 
Transfer, »=1 and « near Unity 
A method of obtaining the enthalpy distribution is first developed, and then 


it is shown how this result may be used to find the approximate position of 
separation in an adverse pressure gradient case. 


130 


BOUNDARY LAYER SEPARATION 


3.1. Enthalpy distribution 


From equations (11) and (12) it follows that with w=1 the enthalpy equation 
(5), with the independent variables changed from x,u to x, y=u/u, becomes 


On u,? On 01? u, dx on 


dx u, On 


Let l=h+ 


+1 ‘ ‘ (5) 


I’—> 0 as o—>1 with no heat transfer. Also, the boundary conditions on J’ with 
no heat transfer are 


or’ 


On = 0, y= 0: y= 


The equation is linearised, i.e. « is assumed to be sufficiently close to 1 for all 


product terms involving (1—«)? and so on, to be ignored. With -, written for 
with «=1, the linearised equation is 


(since in the last term ua? and the I’/I, term is 
ignored because it would appear in a product with 0/’/d). Hence 


ar _ py 


Y 
07, or 


By the Stewartson transformation (14), 7, is defined in terms of the viscous stress 
7, in an incompressible boundary layer with a specified external velocity distribution 
U,(X). If du,/dx is everywhere negative this incompressible case can be solved by 
the approximate Howarth method’. 
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This gives 


[= 
va 
0g PU ,/dxX* 
where ax dU D 


and t(y,) and D(é) (x/(dx/dé) in Howarth’s notation) are known functions, the 
former having been calculated approximately from Howarth’s results (see Tables II 
and III). Hence 


(147 ~ Me I, Papalty” dx 


du, dx ual, /Ta) 


and — 
dx 


(18) 


Accordingly, when du,/dx is constant equation (16) is approximately 
equivalent to 


ar 


(1+ 


(3y—2)u,2M2D al’ 

ag 
(19) 
with the boundary conditions 


or’ , 
On Q, 0; 0, n=l. 


At the leading edge of the boundary layer €=0 and, if u,—u, there, equation 


(19) becomes 
| 


from which /’ at €=0 can be found by numerical integration (cf. Table V—pp. 
138, 139). 


Also equation (19) when differentiated with respect to € at €=0 becomes 


2 (1+ 


M, | 
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TABLE II 
t(, 


€=0 0.0125 0.025 0.0375 0.050 0.0625 0.075 0.0875 0.100 0.1125 0.120 


n=0} 0.332 0.310 0.287 0.263 0.238 0.211 0.182 0.149 0.109 0.062 0 
0.05 | 0.332 0.312 0.292 0.271 0.248 0.225 0.201 0.177 0.150 0.121 0.108 
0.10 | 0.331 0.314 0.296 0.277 0.257 0.237 0.218 0.198 0.178 0.158 0.150 
0.15 | 0.331 0.315 0.299 0.283 0.266 0.248 0.233 0.216 0.199 0.186 0.179 
0.20 | 0.330 0.316 0.301 0.288 0:273 0.257 0.245 0.230 0.217 0.207 0.202 
0.25 | 0.328 0.316 0.302 0.291 0.278 0.265 0.256 0.243 0.231 0.224 0.220 
0.30 | 0.326 0.315 0.302 0.293 0.281 0.271 0.261 0.253 0.242 0.238 0.234 
0.35 | 0.322 0.312 0.301 0.293 0.233 0.276 0.267 0.260 0.251 0.247 0.244 
0.40 | 0.316 0.309 0.299 0.292 0.282 0.277 0.269 0.264 0.258 0.253 0.250 
0.45 | 0.309 0.303 0.295 0.289 0.280 0.276 0.270 0.265 0.260 0.255 0.253 
0.50 | 0.300 0.295 0.288 0.283 0.276 0.272 0.268 0.263 0.258 0.254 0.254 
0.55 | 0.290 0.285 0.280 0.274 0.268 0.265 0.261 0.258 0.254 0.250 0.251 
0.60 | 0.276 0.272 0.268 0.262 0.259 0.255 0.252 0.249 0.246 0.243 0.244 
0.65 | 0.260 0.256 0.253 0.248 0.246 0.241 0.239 0.237 0.235 0.232 0.233 
0.70 | 0.240 0.236 0.234 0.232 0.230 0.226 0.224 0.221 0.220 0.218 0.218 
0.75 | 0.217 0.214 0.212 0.210 0.209 0.207 0.206 0.202 0.201 0.199 0.198 
0.80 | 0.191 0.188 0.186 0.185 0.183 0.183 0.182 0.178 0.178 0.175 0.175 
0.85 | 0.158 0.157 0.157 0.156 0.154 0.152 0.150 0.148 0.147 0.146 0.147 
0.90 | 0.120 0.120 0.120 0.119 0.118 O.115 0.114 O.111 O.110 0.110 
0.95 | 0.071 0.070 0.070 0.070 0.069 0.067 0.067 0.066 0.065 0.064 0.064 
1.00 0 0 0 0 0 0 0 0 0 0 0 


TABLE Ill 

Dé) 
0 0.0125 0.025 0.0375 0.050 0.0625 0.075 0.0875 0.100 0.1125 0.120 
D 0.000 0.024 0.046 0.066 0.084 0.100 0.115 0.128 0.138 0.147 0.151 
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since from equations (18) 


=) = _ 1 (du, 

dé d Uy (1+ 

This equation can be solved numerically for 0/’/0€ at ee (cf. Table IV). 
Elsewhere /’(€, n) can be estimated roughly by making the crude assumption that 


I’ varies parabolically with €, so that 


_ 2 ol’(0, 


Accordingly at € equation (19) will become 

3y 


») 


M.? 
This equation can be solved numerically for J’, so enabling /’ to be estimated for all 
values of € by means of the assumed parabolic distribution. The accuracy of the 
method can be assessed by solving equation (22) at two values of €, say 0-0625 and 
0-120 (the separation value), so obtaining two different estimations for the 
distribution of /’ with €, the maximum discrepancy between which, occurring at 
£=0-12, will be a measure of the maximum error involved in the method. 


In practice the integration of equations (20), (21) and (22) presents some 
difficulty. The standard method of solving such equations, of the form 


* +Qf=R (P,Q,R functions of 7) 


involves finding the solution of 


07h Poh 


ont On +Qh=0, 


but here A cannot easily be found by numerical integration for the range 0-9 <  < 1 


owing to the singularity in ¢ at »=1. However it can be established that at ¢<—0 
and near »=1, ¢ and » are related by an equation of the form 


(1-1) =AV=| 1 (10s *) | 
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TABLE IV 
ar. 
ag AT €=0= 
1+ 7M 
V W 
0 0.504 0.0318 


0.85 0 0.0101 
0.90 - 0.036 0.0045 
0.95 —0.048 —-—0.0011 
1.00 0 0 


and by fitting this relation to the known values of ¢ fairly near to »=1, A can be 
determined. The expression permits analytical integration for /’ and 0J’/d€ near 
n=1 at €=0. For £40, however, such analytical integration is rather more difficult, 
and an easier approximate method of solving equation (22) is as follows. 


The standard method gives a relation for /’(£) of the form 


(€)=L+MI (¢,0-9) 
135 


0.05 | 0.502 0.0318 
0.10 | 0.495 0.0319 
0.15 | 0.483 0.0320 
0.20 | 0470 0.0321 
0.25 | 0.451 0.0323 
0.30 | 0.430 0.0325 
035 | 0.403 0.0326 
0.40 | 0.379 0.0325 
0.45 0.349 0.0321 
0.50 0.313 0.0315 
0.55 0.276 0.0303 
0.60 0.235 0.0286 
0.65 0.192 0.0263 
0.70 0.144 0.0235 
0.75 0.096 0.0200 
0.80 0.045 0.0153 
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where the functions L and M can be found for the range 0 << » <0°9. Also I'(é, 0-9) 
can be found approximately from the condition 


1 


0 


(which can be shown to follow from equations (4) and (5) for any (positive)o with 
© =1 and zero heat transfer), by means of the assumption 


Thus /’ (€), and hence the distribution of /’ with €, can be approximately determined. 


Solutions calculated as described for the cases M,=0, 4, and °© appear in 
Table V. It will be noted that the discrepancy between the two estimated sets of 
results for /’ at €=0-12 is, even in the worst case M,—> 00, quite small, and this 
suggests that the results obtained by the apparently crude method of estimation are 
probably not very much in error. 


3.2. The separation position 
From equations (8) and (15) 


where p’=pal’/I, and —>0 as o—> 1. 
Also with o=1 


P= Papa (7) (equation (12)) 


Hence equation (4) becomes, in terms of independent variables x, », 


Y 
2 ( y-1 )! ar 1 ar 
an? pa (1 — 1+ Me? 3 + papal, ly 


dp, 


with the boundary conditions 


u,*\y-1,,. 
pon [is 
q=z1, r=0; =F 
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Hence integrating, 


" 


n 0 


1 4 
2 
0 


0 


» 0 
A solution for 7 which approximately satisfies this equation and the boundary 
conditions is found as follows. 


It is required to find the separation point, i.e. the position at which 7 (0, x) 
becomes zero. 


Consider a boundary layer where o=1, and where the velocity outside it at 
distance x from the leading edge is equal to the velocity u, outside the boundary 
layer with o near 1 under investigation at distance x from the leading edge. It 
must be possible to choose the external velocity distribution for the «=1 boundary 
layer (i.e. to choose x as a function of x) so that, if 7 (0,x) is the wall stress at 
x in the o=1 case, 

7 (0,x)=T (0, x) 


for vanishingly small x, and so that separation with the ~=1 boundary layer occurs 
at the same external velocity as with the o near 1 boundary layer. Thus for all x, 
z (0, x)=F (x) T (0, x), where F=1, x=0, and F is everywhere finite and >0. If x 
and F could be determined to satisfy this relation, this would give the solution to 
the problem, for a boundary layer for which e=w=1 can be solved, and hence the 
point where T (0, x), and therefore = (0, x), became zero could be found. Accordingly, 
appropriate relations for x ard F are looked for. 


Making the trial assumption 


7 x)=F (x)T x), . (24) 
dr dx aT 
ax dx dx ox’ 


and since the pressure gradient in the o near 1 case is dy/dx times that in the c= 1 
case (because dp/dx= — pa u,du, / dx), equation (23) with its equivalent 
for the «=1 boundary layer gives 


1 4 
+ j F  &) 
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The functions 


n 0 


n 0 


are, when plotted against 7, of roughly similar shape, so that if equation (25) is 
satisfied when »=0, say, or when integrated from 7 =0 to 1, it will be crudely satisfied 
for all ». Thus assumption (24) will be roughly consistent with the equations of 
motion if 


0 


or, since, as can be shown to follow from equations (4) and (5), 


if 
1 


1 
3 I, y-! dF n (1 — 


where 7,, a more convenient approximate equivalent of T, is the viscous stress in 
a boundary layer with o=1 and the same external velocity distribution u,(x) as 
the case under investigation. An F term appearing in each equation as a multiplier 
has been omitted since F=~1 if o is sufficiently close to 1. 


Equations (26a) and (26b) are of the form 


dF 


The values of A and B depend on which equation is used: from the size of the 
(fairly small) discrepancies can be judged how closely assumption (24) is consistent 
with the equations of motion of the two boundary layers, 
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Assumption (24) should also be consistent with their boundary conditions, 
so that 


on 
OT (0, x) _ dp y-1 
2 
and F°*T (0, x) Ba ¢ + +p «| Uu 
Hence dx — FF? = (28) 
dx 


Pa (1 + M.*) + Pa (1 + M.*) 


Unfortunately this is inconsistent with equation (27) at the leading edge x=0, where 
A and B are both proportional to x if the external velocity at x=0 is not zero, 
whereas pu’, does not vanish. Hence some different relation between dy/dx and F 
must be assumed which will satisfy the boundary conditions approximately. One of 
the simplest assumptions would be to take F=1, but although this satisfies the 
boundary conditions at separation where A= — p’y/[ua(1+4 (—1)M,*)], it does not 
agree with them very well elsewhere. Instead, it is assumed that 


n(a) 


n(x)=1- 


where c is a constant, 0<c< 1, and x, is the separation distance of the 
o==1, u,(x) boundary layer. Equation (27) accordingly becomes 


dF A — pw 
aif . (30) 
dx B y-l 

+1 Mi ) 


so that if c is taken equal to 0 and n=1, dF /dx would be 0 and F would be 1. But 
with F=1 the boundary conditions are not satisfied sufficiently closely. Hence c must 
be taken >0 and close to 1 so that relation (29) does not differ greatly from (28) 
except near the leading edge. However c cannot be taken too close to 1 as for 
c=1 relations (28) and (29) become identical and (28) is inapplicable at the leading 
edge. This might seem to imply that the choice of c must be unsatisfactorily 
arbitrary, but this is not entirely the case. For if equation (30) is integrated for 
various values of c to give F (c) at, say xu, it is found that the values of F obtained 
increase fairly steadily until c reaches the region of 0-9, after which they rise sharply 
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to 00 when c=1. This sharp rise is the consequence of attempting to satisfy the 
boundary conditions too closely in the vicinity of the leading edge, and is clearly 
fictitious. Hence it would seem to be a fairly satisfactory compromise to 
choose c=0°9. 


With n(x) thus determined the functions A and B remain to be estimated before 
equation (30) can be integrated to give F. Now if du,/dx is everywhere negative, 
the Howarth-Stewartson method can be applied to the c=w=1, u,(x) boundary 
layer, and gives 


where 


dx u, 


y-1 
dé du, |dx (14 ,| (3y -- 2) 


D(é) 


4 


(equations (18)) 


Hence € can be found as a function of x, and 


i i 
[(é). 
,_ if equation (26a) 
| —»)dn is used 
ui [7] (‘) 
if equation (26d) 
| (1— 1?) dn is used 
t 


where 


Also, since the boundary conditions give( 1 
N/w 


Ot /On/] w 
= if (6a) is used. 


1+ 


{ tar if (26b) is used. 


2 
(0t/On)w 
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TABLE VI 


i ; 
AVERAGE VALUES OF AND -dn AGAINST & 
w 


t 


0 0.025 0.050 0.075 0.100 0.1125 0.120 


0.754 0.844 0.983 1.228 1.93 3.26 0 


0.265 0.275 0.293 0.322 0.378 0.425 0.467 


i 
At = 1.63 


Strictly speaking i, found as described earlier, is, for « near 1, a function of », €,M, 
and the external velocity distribution, but for the purposes of the present method 
i can be taken as a function of » and € only, and thus in Table VI approximate 


average values of 
0 


i i 
1 


are presented as functions of €. The integrals 


1 1 
10-9) nd | 

0 


0 


are given in Table VII. 


can be readily found from Table II and the condition (t dt/@n),=£, and finally 


1 


2 
fray 


0 


is given in Table VII. Hence for any Mach number and any external velocity 
distribution with du, /dx everywhere negative, F and dx /dx can be found quite easily 
from equations (30) and (29). 


Having calculated 0x/0x with dy/dx determined, the change of separation 
distance as compared with the == 1, u, (x) boundary layer (which will have been 
solved by the Howarth-Stewartson method) could be found by solving in the same 
way the e=w=1, u, (x) case corresponding to, say, s=0-9. However, the changes 
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TABLE VII 


1 


1 1 
n (1-1) j n (1-7?) 2 
dn, dy AND (@t/anw t dn AGAINST & 
0 0 


0 


£ 0 0.125 0.025 0.050 0.075 0.100 0.1125 0.120 


ae 0.664 0.680 0.696 0.725 0.765 0.804 0.823 0.838 


[22-6 1.028 1.047 1.063 1.112 1.162 1.213 1.243 1.259 


0 
1 


co 612.59 5.65 2.19 1043 0.434 0.212 0 


involved are small and it would require very accurate computation to find them at 
all precisely. A better method of determining them is as follows. 


According to an approximate method due to Thwaites‘, separation for an 
incompressible boundary layer with external velocity distribution u,(x) occurs 
at x, where 


Hence for the compressible boundary layer with «=1 and external velocity 
distribution u, (x), separation occurs where 


-0.182( 


x 


= 


0 


— 0.182 


where by the Stewartson transformation 


hy" 
U,=u, 
d dX _ 1) 
an an Ia 


(cf. equations (14)). Thus at separation 


2 


du, Ri 
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1 


(1, 


If =C (x) 


du,/dx 1+ yl 
2 
C C (Xm) + (Xs — Xn) (AC 


where xs, is the o=1 separation distance, and the change of separation distance 
Xs— Xs, is assumed to be small. Also 


—0.182 C dx 
d 


(since as e—> 1, dy/dx—> 1), and 


dx dx dx dx 
(dx /dx)z, = (dx/ dx)z,, 


(32) 
where P= 


As a test of accuracy the s=0, »=1 case, for which the correct solution can 
be found as described earlier, can also be solved by the present method if the Mach 
number is small. For then »’ is small and F, dy/dx~1, T ~7, as assumed. 


With «=0, 


B | 
2 
d 


and the terms of equation (30) can be calculated readily if equation (265) is used. 
The change of separation distance found in this way for the case with u,=ua— x, 
8 constant, is 75 per cent. of the true change, for the true separation point values of 
u,/Ua With c=o=1 and c=0, w=1 are found to be (0°880+0-0092 M,”) and 
(0-880 —0-0066 M,”) respectively, while the value calculated by the present method 
for the «=0 case is (0:880—0-0027 M,7). Thus the method works satisfactorily 
in this instance. 
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w=1, o near | are, at M,=4, 


Xs— 


0-28(1—«) 

and at M,=10, 
_ 9.96 (1-0) 


The results given by the method for linear external velocity gradient cases with 


if equation (26a) is used. 
if equation (26b) is used. 


if equation (26a) is used. 
if equation (26b) is used. 


It will be seen that little difference is made to the predicted change of separation 
distance if equation (265) is used instead of (26a). The differences in dy/dx and F 
are correspondingly small, and this means that the assumed solution fits the 
equations of motion fairly well. The boundary conditions are not seriously violated 
either, for the average value of 


Pha (1 
Q 


(dy | dx — F?) 


over the whole boundary layer is in the region of 0-9 and O0<Q<1: if QO were 
everywhere equal to 1 the boundary conditions would be satisfied accurately (apart 
from errors of order (1—«)?). 


Hence it seems safe to conclude that when o is near 1 the proportionate 
increase of separation distance, compared with the o=1 case, is small at moderate 
Mach numbers (0-5), and probably not much greater than (1—«) even for Mach 
numbers of 10. 


This conclusion receives support from Young’s approximate method"? based 
on the momentum equation. Young obtains the following general relation for 
separation distance x, 


r 


0 
a ot| 


—3 PU," pa 
f (du, / dx) py 


f=9-072 [1+ 0-365 (y— 1) M,”]'-° 


g=9'18+ 1-436 


1 
Bia 
Pi=Pa 1, 


and it is assumed that 


[14 > M, 1) 


Differentiation of this relation with respect to o at c=w=1 for the linear external 


(33) 
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i velocity gradient case gives that, with Mach numbers of 4 and 10, the proportionate 
increases in separation distance are 0:28 (1—«) and 0-41 (1—c) respectively. 


4. Approximate Solution for Boundary Layers with no Heat 
Transfer, c=1 and » near Unity 


With o=1 and no heat transfer the equation of motion plus continuity (4) 
becomes, in terms of variables x, y=u/u,, 


| = +[o (y—1)+y] (equation (13)) 


with the boundary conditions 


If | 1—©]| is sufficiently small 


+(1 | may be replaced by + (1 o)log| 1 +(1 | 


and hence 


1+(1—wo)log| 1+ M, } 
Or 


I, 
=—— (142 SiMe Ja 2 Oy + papa,” 


[o(y— 1) +y] Ma Ma 


dx 2 On I, ‘ 0 
0 
1 


id 


=a 
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Consider the trial solution 


7 (yn, x)=F (x) T (n, x) 


where T is the viscous stress in a sc=w=1 boundary layer such that the velocity 
outside it at x is equal to the velocity u, outside the w near 1 boundary layer at x. 
By arguments similar to those used for the »=1, o near | case, it can be shown 
that the boundary conditions and equation (34) will be approximately satisfied if 


F=1, x=0, 


—--—--, §,=S at x=0, 


wn 
ill 
| 
| 
| 
| 
| 


(7, +07, an). 


1 
(y— 1) pau,’ dy 


1 

log 

R= n=1-(2) , c=0°9, 


7, is the stress in the w=o=1, 
distance is x,, and 


tog (1+ mr) 
ap S+R-S,) -1 | (7, +07, 


u,(x) boundary layer for which the separation 


S+R+S, 


Ty 


(The relation for dx/dx—F* would satisfy the boundary conditions accurately, 
apart from errors of order (1—)*, if m were zero). 


According to this solution F is everywhere finite and >0 so that the value of 
u, for which T (0, x) is zero is also the value for which 7(0,x) is zero, i.e. the 
separation point. This can be determined by solving for dy/dx and using equations 
(31), (32) (which, although developed for the w=1, o near 1 case are equally 
applicable here). 
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Now if du,/dx is everywhere negative 


du,/dx 


dé_ du, | (equation (18)) 
where —= +} 
dx 
(3y —2) M,?u,du, / dx 
ua(I, / Ta) ] 


so that R=QV where 1 
M2? - 
Q= 
1+ yt (t+ 0t/On)w 
(7, +07,/On)w (1+ Ma ) du, 
= 
1 
and S=log (1+ mt) —-W(Q,&), 
since log[1+ (i |- log (1+ me) + log (1 7°Q). 
TABLE VIII 


GeGi G2 GS @ 7 08 09 095 1.0 


€=0 | 0.040 0.083 0.129 0.180 0.233 0.292 0.358 0.430 0.518 0.570 0.671 
0.025 | 0.036 0.075 0.120 0.165 0.215 0.272 0.333 0.400 0.480 0.527 0.615 
0.050 | 0.031 0.065 0.101 0.143 0.187 0.235 0.289 0.348 0.416 0.458 0.530 
0.075 | 0.025 0.052 0.081 0.114 0.148 0.185 0.227 0.272 0.324 0.357 0.410 
0.100 | 0.016 0.034 0.051 0.071 0.092 0.115 0.139 0.166 0.196 0.215 0.246 
0.110 | 0.011 0.020 0.035 0.048 0.060 0.076 0.093 0.111 0.131 0.143 0.151 


W (Q, 0.120) =0 W (0,£)=0 
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V (€) can be calculated from Tables VII and II and the relation (¢ 0t/0n)~=£€, and 
W (Q, €) is given in Table VIII. Accordingly, for any Mach number and any 
external velocity distribution with du, /dx everywhere negative, dx /dx, and hence the 
change of separation distance, can quite easily be calculated once € has been found 
as a function of x from equations (18). 


It was thus calculated that with Mach numbers of 4 and 10 and constant 
adverse external velocity gradients, the proportionate increases of separation 
distance compared with when »—1, are 0-18 (1—w) and 0:53 (1—w) respectively. 
The corresponding results given by Young’s method (Ref. i, cf. equations (33)) 
are 0:17 (1—©) and 0:26 (1— ©). 


5. General Conclusion 


From the conclusions reached in Sections 3 and 4 it may be inferred that, with 
the types of pressure distribution most usual in practice and free stream Mach 
numbers up to 10, the laminar separation point in compressible flow with zero heat 
transfer and 0°7<0, » <1 probably occurs at roughly the same position as with 
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Buckling of Oblique Plates with Clamped 
Edges under Uniform Compression 


W. H. WITTRICK, M.A., Ph.D., A.F.R.Ae.S. 


(Department of Aeronautical Engineering, University of Sydney) 


Summary: The problem of the buckling of flat oblique plates (i.e. in the shape 
of parallelograms) with all edges clamped under a uniform system of stress is 
considered. An energy method of calculating the critical stresses is developed and 
applied to the particular case of a clamped oblique plate under uniform compression 
parallel to two of the sides. Numerical values for the buckling stress coefficient are 
obtained for plates having acute angles of 60° and 45°, and side ratios varying 
between 3/5 and 5/3. 


1. Introduction 


In the past, the skin panels of aircraft have been mainly rectangular, or very 
nearly rectangular, in shape. The advent of swept wing aircraft, however, has 
introduced the possibility of oblique panels, which are bounded by parallelograms. 
A certain amount of work has been done on the behaviour of such panels under 
normal loading but, apart from the work of Anderson’, Guest” and that described 
in the present paper, no results appear to be available for the stability of oblique 
panels. 


Anderson considered the problem of a flat sheet, extending to infinity in all 
directions, subdivided by a series of equally-spaced non-deflecting supports into an 
array of identical oblique panels, and subjected to compressive stresses both parallel 
and perpendicular to one support direction. In the corresponding problem for an 
array of rectangular panels the nodal lines of the buckled plate are straight and 
parallel to the sides of the panels, so that adjacent panels buckle identically but in 
opposite directions. The boundary conditions for an individual panel, therefore, 
are those of simply-supported edges. In the oblique case, however, the nodal lines 
are in general curved so that interaction between adjacent panels must occur. The 
problem is then similar to that of an array of rectangular panels in shear“ and 
cannot be treated in terms of an individual panel. 


Anderson simplified the problem considerably by assuming that the buckle 
pattern repeats itself over certain numbers of panels in both directions, and that any 
nodal lines which may occur are straight, although not necessarily parallel to the 
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sides of the panels. The expression which he adopted for the mode contained four 
unknown parameters which were adjusted to give the lowest buckling stresses by 
the Rayleigh-Ritz method. In a few cases he used the basic features of this 
approximate mode to check the accuracy by a more accurate analysis and obtained 
quite good agreement between the two. 


Anderson’s results were presented in the form of graphs of buckling stress 
coefficient against aspect ratio for a series of angles of sweep, and for each of the 
two compressive stresses acting alone. He also gives interaction curves for equal 
sided panels under the combined action of the two stresses. 


Guest™ considered the problem of an oblique plate with all edges clamped 
under uniform compression parallel to two of the sides by the method of Lagrangian 
multipliers which gives both upper and lower bounds for the buckling stresses. 
Numerical results were obtained for the lower bounds for two equal-sided panels 
with angles of sweep of 30° and 45°. 


The present paper considers the problem of an oblique plate with all edges 
clamped under any uniform system of edge stresses. The Rayleigh-Ritz method 
is used and the problem is reduced to finding the roots of determinants which are 
upper bounds to the true buckling stresses. 


The method was applied numerically to the particular case of clamped oblique 
plates under uniform compression parallel to two of the sides. Values of the side 
ratio varying between 3/5 and 5/3, and two different values for the acute angle of 
the parallelogram, namely 60° and 45°, were considered. The results for the 
buckling stress coefficient are given in Table I and plotted in Fig. 4. In order to 
obtain these results it was necessary to evaluate approximately fifty eighth-order 
determinants. The results were in good agreement with the two lower bounds 
determined by Guest. 


It is, of course, realised that clamped edges are not the most practical, but they 
are by far the easiest to deal with. It is for this reason that attention has so far 
been confined to them. 


Notation 


a,b lengths of the sides of the parallelogram 
=a/b 


Se 


angle of “sweep” of the parallelogram, its internal angles 
being +) 


thickness of the plate 
Young’s modulus 
Poisson’s ratio 
flexural rigidity of the plate = Er’ /[12 (1 — 
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x,y Cartesian co-ordinates (Fig. 1) 


u,v _ oblique co-ordinates (Fig.1) 


€,» non-dimensional oblique co-ordinates equal to u/a, v/b 
respectively 


w lateral deflection of the buckled plate 
fe» fy, fey components of the applied stress system (Fig. 2) 
B. =f,b?t costa/D 
Bay =f tcos'a/D 
B, =f,b?tcos?«/D 
A,,4, determinants, one of which is zero when the stress system is 
critical 
Ampnq the general element of the determinants A,, A, 
Q,, coefficient of the general term in the series representation of 
w (equation (10) ) 
F,,G, functions defined by equations (12) 
definite integrals defined by equations (16) 
m,n,p,q,P,Q,r positive integers 


2. Basic Theory 


Consider a thin flat plate, of uniform thickness t, bounded by a parallelogram, 
as shown in Fig. 1. Let the plate be subjected to a uniform stress system f,, f, 
and f,,, as shown in Fig. 2, and suppose that the magnitudes of these stresses are 
such that buckling has occurred. The plate is then in a state of neutral equilibirum 
and the total potential energy of the system is a minimum. If every edge of the 
plate is either clamped or simply supported this condition may be expressed by 
the equation” 


WS +55) +2 )|awdxdy=0. (1) 


where D is the flexural rigidity of the plate, w the lateral deflection at the point 
whose Cartesian co-ordinates are (x, y), and 45w is any variation in w which satisfies 
the boundary conditions. 


In order to simplify the equations of the boundaries it is convenient to transform 
equation (1) from Cartesian co-ordinates (x, y) to oblique co-ordinates (u, v) as 
shown in Fig. 1. The equations governing the transformation are 


u=x—ytana, v=yseca ‘ P (2) 
and the boundaries of the plate are now given by u=0, a and v=0, b. 
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Vv 
me) 
| 
qd 
Fig. 1. 
Dimensions of the plate. 
Further we have 
0 0 0 


Equation (1) then becomes 


ab 
00 
2 2 
2 
+f, cos? —2sin a +sin?a ] dw dudv = ‘ (4) 


Finally, on writing 


the boundaries of the plate are given by €=0, 1 and »=0, 1, and equation (4) 
becomes 


11 
4 4 4 4 4 


+ +2Bay(2 + 
and B,=f, b*t cos*a/D 
Bey = fryb*t cos*a/D 


8B, =f, cos?«/D 
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Now, since 4w is any variation in w which satisfies the boundary conditions, it 
follows that the expression in square brackets in equation (6) must itself be zero. 
This gives the differential equation for the deflection w, but it cannot be solved 
exactly. In what follows only the case of a plate with all edges clamped will be 
considered, and approximate solutions will be obtained by the Rayleigh-Ritz method. 


3. Oblique Plate with All Edges Clamped 


So far the equations developed are applicable to a plate in which every edge 
is either clamped or simply supported. The simplest case to deal with numerically 
is when all the edges are clamped, and the analysis is now restricted to this 
special case. 


The boundary conditions of zero deflection and zero normal slope at every 
edge can be written 


ow 

on =0 at n=0, 1 

(9) 
at é=0, 1 


To obtain numerical values for the buckling stresses, w is now expressed by 
the series 


iy 


| 
“Tey 


Fig. 2. 


Definition of applied stress system. 
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where, in order to satisfy (9), 


= 1 = => 
F, (0)=F, (1)= F,(0)=F,(1)=0 


G, (0)=G, 


and the dashes denote differentiation. 


Several possibilities exist for the choice of the F and G functions, but in the 
numerical work described in this paper the functions used are 


F, ()=€ +(- sin pre 
(12) 


Z 
G,(n)=n +(- 7? qm 


The functions (12) were first used by Iguchi’ for the determination of the buckling 
load of a clamped rectangular isotropic plate, and later by Smith in the 
corresponding problem for an orthotropic plate. The boundary conditions in both 
these cases are identical with (9). Smith also discusses the validity of the 
expansion (10). 


If equation (10) is substituted into (6) and the choice is made that 


bw =F, (€)G, (n) ‘ 
the set of equations 


p=1 q=1 2, 3,....Q. 
is obtained, where 
Awana = +X? Lng +2 (1-+2 — 
4 sin +2 Lng ) + Be + 
+2 Bey Ing’? —Sin Ing nq!) + 
+ By (2 Sit Imp + (15) 
and Ing = | Fn Fr OdE= | 16) 
: 


It is easy to see from equations (12) that F,(€) is symmetrical about =} if 
p is odd, and anti-symmetrical if p is even, so that 
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It then follows, from the definition (16), that 


if (m+ p) is odd and r even, or if (m-+p) is even and r odd. It is seen therefore, 
from equation (15), that Amp =O if (m+p+n+q) is odd. Consequently any 
equation in the set (14) in which (m+n) is even contains only those coefficients a,, 
in which (p+q) is even, while any equation in which (m+n) is odd contains only 
those coefficients a,, in which (p+q) is odd. The set of equations (14) therefore 
splits up into two distinct sets. The first set contains only the a@,, in which (p+ q) is 
even, while the second set contains the a,, in which (p+q) is odd. Now all the 
equations in either set are homogeneous in the a,,, so that a non-zero solution can 
exist only if the determinant formed by the coefficients of the a,, is zero. Therefore 
the following alternatives exist : — 


Either (i) the plate buckles in a mode in which all the terms of the series (10) 
correspond to (p+q) even, and the determinant 4, formed by the 
coefficients of the corresponding a,, in the set of equations (14) is 
Zero, i.e. 


m,n 


where both (p+q) and (m+n) are even, 


or (ii) the plate buckles in a mode in which all the terms of the series (10) 
correspond to (p+q) odd, and the determinant A, formed by the 
coefficients of the corresponding a,, in the set of equations (14) is 
zero, i.e. 


where both (p+ qg) and (m+n) are odd. 


Since the elements Ampn, of the determinants A, and A, are functions of 8,, 8, and 
By, equations (19) and (20) are equations from which the buckling stresses may 
be determined for any particular values of A and z. The question of which of the 
two equations is the one to choose depends on which one gives the lowest buckling 
stresses and can only be resolved numerically. It will be seen later that for some 
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Plate under uniform compression. 


values of A and 2 equation (19) gives the lowest stresses, while for other values 
equation (20) is the one to choose. 


Considering the buckling modes in more detail, it follows from equation 
(17) that 


F, (€)Gq(n)=(— 1? ** F, (1—£)G, (1-1) (21) 


Hence, from equation (10) 


w(&,n)= w(l1—€,1-n) if (p+q) is even . ‘ (22) 
= —w(l1—£,1—n) if (p+ q) is odd (23) 


It is to be expected that the buckling mode will in general contain a number of 
nodal lines, and it is now clear from equations (22) and (23) that if the plate buckles 
under stresses given by equation (19) the mode will contain an even number 
(including zero) of nodal lires, while if the buckling stresses are given by equation 
(20) there will be an odd number of nodal lines. Just how many such nodal lines 
there will be in any given case is impossible to determine without resorting to 
computation. 


Before considering specific cases it is noted that, by virtue of equations (15) 
and (18), if (m+ p) and (m+ q) are even, 


1 
| +A? + 2 (1 + 2 sin?a) Ting + 


= 


+ Be Imp"? Ing? — 2 Bay Si & Img"? Ing’ + By (A? Ing’? + ing?’ Ing”) 
(24) 


while if + p) and (n+ q) are odd 


Angug= + ) +2 ag! By Tag! Mag sin 2 
(25) 
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Also, on using equations (11), and integrating equation (16) by parts, it can be 
shown that 


From this it follows that 


A mpnq — A pmqn (27) 


which means that each of the determinants A, and A, is symmetrical about its 
principal diagonal. 


Expressions for the various definite integrals /,,,“? which occur on the right 
hand sides of equations (24) and (25) are given in the Appendix. 


4. Clamped Oblique Plate Under Uniform Compression 


The special case is now considered of an oblique plate with all edges clamped 
under uniform compression parallel to the sides of length a, as shown in Fig. 3. 


Then ; ; (28) 
and equations (24) and (25) give 
(29) 
if (m+ p) and (n+ q) are even, and 


if (m+ p) and (n+q) are odd. 


Using equations (29) and (30), and the-expressions for the various definite 
integrals given, the determinants A, and A, were constructed. For this purpose the 
values of P and Q, defining the number of terms in the series (10), were both taken 
to be equal to 4, so that the determinants A, and A, were both of eighth order. This 
gave what will be referred to as the master determinants, in which the elements con- 
tained A, « and 8, in symbolic form. Two values for the angle 2 were considered, 
namely 30° and 45°. These, and several different values of the side ratio A, varying 
between 3/5 and 5/3, were inserted in the master determinants, and in each case the 
lowest value of 8, required to make the determinant zero was calculated The method 
used was as follows. First a value of 8, was estimated, using as a basis the known 
results for a rectangular plate («=0) and any previously computed results for the 
oblique plates. This value was inserted in the determinant which was then evaluated 
using Crout’s’ method. In all cases it was found that if the determinant was positive 
the estimate of 8, was too small, while if the determinant was negative the estimate 
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was too large. Accordingly a second estimate of 8, was inserted in the determinant, 
which was again evaluated. Using these two results a third estimate was made, the 
determinant being evaluated once more. Finally, from the three results a parabolic 
interpolation fixed the required root to at least four-figure accuracy. The critical 
values of (8,/*) found in this way are given in Table I and plotted in Fig. 4, which 
also shows the results for a rectangular plate («=0) as determined by Levy“. 


TABLE I 


CRITICAL VALUES OF (8,/77) 


r 3/5 4/5 1 5/4 3/2 5/3 


From A,=0 12.67 8.81 7.67 7.41 


From A,=0 9.98 6.44 5.41 5.28 -— >i 


From A,=0| >9.98* — 617 5.14 485 4.90 


The two results marked by an asterisk in Table I were obtained in the 
following way. In order to show, for example, that for 2=45° and A=3/5, the 
equation 4,=0 gives a higher value for 8, than A,=0, the value 8,=9.98"? was 
substituted into 4,, which was then evaluated and found to be positive. This fact 
means that the lowest ‘root of A, is greater than 9.987. 


It is known that rectangular plates with A=3/5 and 5/3, buckle in modes 
corresponding to (p+q) even and odd respectively. It has been shown that the 
same applies to oblique plates with «=45°, and it is therefore assumed to be true 
for «=30°. 


Since only the smallest values of 8, are required, the curves drawn in full in 
Fig. 4 give the required critical values. It is by no means obvious which is the 
best non-dimensional form for plotting. It is clear, however, that as A—> oO, 
whatever the angle a, the buckling stress must become equal to that of an infinite 
strip with clamped edges: and the same width perpendicular to the direction of 
compression as the oblique plate. From the first of equations (8), 


_ 
~ costa 


fe 


This may be rewritten in the form 


(b cos «)*t 


(31) 


sec? 


where (bcos) is the perpendicular width of the plate. Hence, as A—> 00, 
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Fig. 4 


since this represents the limit for an infinite strip with clamped edges. The limits 
represented by equation (32) are shown in Fig. 4. 


Some evidence exists as to the accuracy of the results given in Table I. 
Guest‘*’ has obtained lower bounds to the buckling stress when 2=30° and A=1, 
and when z= 45° and A=1, by a method using Lagrangian multipliers. The results 
of the present paper, being obtained by a Rayleigh-Ritz method, are upper bounds, 
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so that a comparison between these and Guest’s results gives some idea of the 
errors. The values of 8,/7? obtained by the two methods are shown in Table II 
which also gives the percentage deviation from the mean of the two results. 


TABLE Il 


VALUES OF 


Lower bound Upper bound Deviation 
(Guest) (Present paper) from mean 


z=30°, A=1 | y 7.67 +0.4% 


a=45°, 5.41 +2.2% 


From this comparison it appears that the results for «=30° are more accurate than 
for <=45°. Also it is likely that the error will increase as A increases so that the 
most inaccurate result is probably that for 2=45°, A=5/3. However, this is 
unlikely to be more than about 4 per cent. in error, judging by the deviations shown 
in Table II. 
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APPENDIX 
Values of integrals involved in the clamped oblique 
plate problem 


The following are expressions for the various definite integrals occurring on 
the right hand sides of equations (24) and (25). 


- 
m* 


+5) + me if m and p are odd, 


a)+3 mint 8 if m and p are even, 


0 otherwise. 
5 if m and p are odd, 


if m and p are even, 


0 otherwise, 
2 
3 if m and p are odd, 


2 
= 0 otherwise. 


1 
=} + if m is odd and p even, 


30 
1 
= 


mt if m is even and p odd, 


= 0 otherwise. 
if m is odd and p even, 


if m is even and p odd, 
= 0 otherwise. 
In the above expressions the term 8 is defined by 
s=1 if m=p 
8=0 if 
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A Note on the Theory of the Constant-Area 
Mixing of Compressible Flows as Applied 
to High-Speed Wind Tunnel Design 


F. G. IRVING, M.Eng., D.L.C. 


Summary: The one-dimensional equations governing the mixing of compres- 
sible flows in a channel of constant cross-sectional area are derived and are used 
to investigate the running of an induction-type high-speed wind tunnel. While in 
practice it is usually difficult to achieve working section Mach numbers greater than 
about 1-8 to 2:0 in this type of tunnel, this simple theory gives no obvious reason 
for this and indicates that higher Mach numbers should be attainable. 


It is also shown that an increase in efficiency is predicted, particularly at the 
higher supersonic Mach numbers, if the injected stream enters the mixing section 
at the same Mach number as that of the working section main stream. 


1. Introduction 


If two streams of a compressible fluid, under given initial conditions, are 


allowed to mix in a channel of constant cross-sectional area, then the final condition 
of the mixed stream may be found by applying the one-dimensional equations of 
momentum, energy and continuity. To simplify the computation, the extra assump- 
tion that the pressures of the entering streams are equal is often made (e.g. Ref. 1), 
but it is not necessary to do so. In Ref. 2 equations were developed similar to those 
given in the present paper, assuming that the injected stream entered the mixing 
section at sonic velocity. In the first part of this analysis this assumption is not 
made, and more general equations are derived. Two conditions are then examined 
and compared. 


(i) When the injected stream enters the mixing section at sonic velocity. 

(ii) When the injected stream enters the mixing section at the Mach number 
of the induced stream upstream of the mixing section (i.e. at the working 
section Mach number in the case of an induction wind tunnel). 


Notation 
density 


velocity 
pressure 
cross-sectional area 


Paper received January 1952. 
[The Aeronautical Quarterly, Volume IV, February 1953] 


164 


p 
A 


MIXING OF COMPRESSIBLE FLOWS 


speed of sound 

Mach number 

temperature (°K) 

isentropic efficiency of diffuser 

ratio of specific heats 

injector slot area ratio= A,/A, 

defined by equations (7) 

mass ratio=(rate of mass flow of induced air)/(rate of mass flow 
of injected air) 

“ideal” mass ratio 

mass ratio corresponding to the finite slot area ratio, 7 


refers to stagnation conditions at entry to working section. 
working section conditions. 


conditions at the end of the working section in the plane of the 
injector slots. 


conditions of injected stream at exit from injector slots. 
stagnation conditions of injected stream at entry to injector. 
conditions of mixed stream at end of mixing section. 
conditions at diffuser exit. 


2. One-dimensional Equations for Constant-area Mixing 
Consider an induction-type high-speed wind tunnel, as sketched in Fig. 1. 
The following assumptions are made:— 
(a) that the flow between “0” and “1” is isentropic, 
(b) that the flow between “4” and “3” is isentropic, 
(c) that at “5,” conditions are uniform across the channel, i.e. mixing is 
complete, 
(d) that A,=A,+A, and A,=A,, 
(e) that the kinetic energy of the stream at “6” is negligible, 
(f) that the ratio of the specific heats is the same for the injected and induced 
streams. 


The equations of continuity, momentum and energy may now be written for 
the mixing process between “2,” “3” and “5.” 


Continuity : 
. 
oz A,/ As. 


P2+ +o (Ds + (1+ 7) + P5757) 
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Fig. 1. 
Schematic layout of induction-type wind tunnel. 
Energy: 


piv, 2 Es + Opsv =(1+ {4 Bs . 


Now there may or may not be a normal shock between “1” and “2.” If 
not, conditions at these sections will be identical. Otherwise it remains true that 


= (4) 
Py + =P2+ 


These relationships may be substituted in equations (1), (2) and (3), together 
with the usual expressions for isentropic flow between sections “0” and “1” and 
sections “4” and “3.” Simultaneous solution of the final three equations gives 


and equation (2) may be written 


My 


y-1 
(1+ M 


where 
=L 
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Hence from equations (5) and (6), M, and p, may be obtained, knowing the initial 
conditions of the injected and induced streams. 


If T,=T,, as is approximately the case for induction wind tunnels, then a,=a, 
and p,/p,=p,/ p>, and equations (5) and (6) become 


M.*) 


Ps Ds 


To define the conditions under which an induction-type wind tunnel will 
operate, the equations (8) and (9) must be solved simultaneously with the equation 
describing the diffuser behaviour. This may be written 


me] 


The condition that the tunnel shall operate is, in general, that 


3. Discussion of Equations (8) and (9) 


If o is very small, equation (9) may be written 


Ps 2 
—(1+yM,”)=L,+L,o Pe) ‘ . 
Ds Pe 
In these circumstances both p,/p, and M, become functions of M,, M, and the 
product «(p,/p,) only, and it becomes unnecessary to consider o and p,/p, 
independently. 


The performance of the tunnel obtained by the replacement of equation (9) by 
equation (12) may be considered as an “ideal” performance. It will be shown later 
that, for given initial values of M, and M, and a given diffuser efficiency 1,, 7 (p,/p.) 
increases as o increases. Other things being equal, the “mass ratio,” defined as 
(rate of mass flow of induced air)/(rate of mass flow of injected air) is inversely 
proportional to o«(p,/p,) and hence it becomes less as o increases. The “ideal” 
mass ratio, corresponding to an infinitely small value of o, therefore represents a 
maximum value (for given values of M, and M,) which is unattainable in practice. 
These “ ideal” values, however, are useful for purposes of comparison and practical 
values corresponding to finite values of « may be deduced from them, as will be 
shown later. 
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In practice, the smallest injector slot area (i.e. the smallest value of «) which 
can be used is limited by the maximum value of p,/p, which can conveniently be 
applied to the tunnel. It is also important to note that, while minimum slot area 
leads to the best performance of the tunnel alone, it will not, in general, give the 
best performance (i.e. longest running time) for an intermittent induction tunnel 
operated from stored compressed air. In the limit, as the pressure p, approaches 
the storage pressure, the running time will tend to zero. For such a tunnel there is 
an optimum value of p,/p, and a corresponding value of « for each value of M,, 
other things being equal, which may easily be found® for any plant of known 
characteristics. 


For the moment it will be convenient to consider the “ideal” performance, 
if need be modifying the figures so obtained at a later stage to allow for finite 
values of co. 


Equation (8) is seen to be a quadratic in M,’. If its roots are denoted by M,° 
and M,”’, then they will satisfy the equation 


(1+yM,* (1+ 


which is the relationship between the Mach numbers on opposite sides of a normal 
shock. The two roots of equation (8) therefore represent supersonic and subsonic 
conditions after mixing. In practice, the subsonic condition is found to obtain and 
therefore only this root will be considered, it being understood that », is the efficiency 
of the subsequent subsonic diffusion process. 


4. Sonic Injection 


For mechanical convenience, it is usual to make the injector slot of an induction 
tunnel in the form of a convergent channel. In such a case, sonic conditions prevail 
at the exit from the slot and then 


and 


This assumes that 


as is usually the case in practice. Equations (8), (12) and (10) are then to be solved 
simultaneously, having inserted these values of N, and L,. A convenient procedure 
is as follows :— 
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for O<M,?<1. 


On the same diagram plot p,/p, for the required values of 1, 
(equation (10)). 


For each value of M, take a suitable series of values of «(p,/p,) and 
hence obtain the corresponding values of 


(eT 


From curve (a) find the s values of M.”. 


For each value of M, and each value of o(p,/p,), together with the 
corresponding M,”, find p,/p, from equation (12). 


Plot p,/p; against M,” on the same diagram as the diffuser curves (b). 
The intersections of these two sets of curves represent tunnel running 
conditions. 


Hence, for each value of M, and 1, read off the corresponding value of 
(p,/P»). 


A portion of this diagram is shown in Fig. 2. 
Now the mass ratio is given by 


(13) 


And hence for the chosen values of 1, the “ideal” mass ratio may be plotted against 
the working section Mach number M, as in Fig. 3. 


For any chosen value of M, and », it is then possible to plot the mass ratio 
corresponding to finite injector slot area ratios, m., against « by proceeding as 
follows (see Fig. 4):— 


(i) 
(ii) 


(iii) 
(iv) 


Over a suitable range of values, plot p,/p,; against M,”. 


Choose a suitable series of values of « (p,/p,) and hence find the corres- 
ponding values of M,? as in (c) and (d). 


Hence, from equation (9) find the corresponding values of (1+¢) p,/ py. 
For each value of o(p,/p,) choose a suitable series of values of « and 
hence from (iii) find the corresponding values of p,/p,. 
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Fig. 2. 


Diffuser and mixing process curves. The figures beside the plotted points are the appropriate 
values of o (p,/p,). In all cases points are plotted at increments of 0:02. The M,=2°4 line 
is omitted for clarity, 
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Variation of mass ratio with slot area ratio for sonic and supersonic injection 
at M, = 1, =0°7. 


(v) Plot, on the same diagram as (i), p,/p; against M,” for the chosen values 
of and (p,/ py). 


(vi) From the intersections of these curves at the chosen values of o with the 
diffuser curve, find the values of M,? which represent tunnel running 
conditions. 


(vii) From equation (13) find the corresponding values of m, and plot m, 
against 7. This curve for M,=1-6 and 7,=0:7 is shown in Fig. 5. 


5. Supersonic Injection 

While it is usual in actual induction tunnels for the injector slots to be con- 
vergent, so that the injected air emerges at the local speed of sound, the pressure 
ratio across the slots is usually considerably greater than the critical. Hence if the 
injector slot were made in the form of a convergent-divergent nozzle, the injected 
stream would be supersonic as it entered the mixing section. For convenience, 
suppose that the form of the injector slot is adjusted so that the injected stream 
leaves the slot at the same Mach number as that of the main stream in the working 
section. Then M,=M,, M, and M, being greater than unity. 
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Then N,=N, and L,=L,, so that equation (8) becomes 


Gi) 


| 


(14) 


Hence M, and M, are related by the normal shock equation, so that, knowing 
M,, M, may be found from the appropriate tables. (The supersonic condition 
M,;=M, is not considered, as before). 

Equation (9) becomes 


The diffuser equation (10) remains unaltered, so solving equations (10) and (15) 
and inserting the condition p,=p,, the running condition is defined by the equation 
Ps 
1 + (2) 2 
lio y-1 
L, (1 + 1 M; 


Hence for chosen values of M, and 1», the corresponding values of 
[1+o(p,/p,)][1+o]~' may be obtained from equations (14) and (16) directly. Ifo - 
is assumed very small, the corresponding values of [o(p,/p,)]~' are the “ideal” 
mass ratios for supersonic injection. The “ideal” mass ratio curves are shown in 
Fig. 6. The mass ratios corresponding to finite slot area ratios may also be found 
directly from the foregoing equations. The values under these conditions for 
M,=1-6 and »,=0-7 are shown in Fig. 5. If-the “ideal” mass ratio is denoted by 
m,, then from equation (16) 


or m= i+o(l+m,) ( l 7) 


where m, is the mass ratio corresponding to the finite slot area ratio 7. Fig. 7 
shows curves of m, against o for various values of m,. 


Consider now the curves of mass ratio against slot area ratio for sonic and 
supersonic injection. The ratio of the mass ratios (which is also the ratio of the 
injector mass flows, since the working section mass flow is the same in each case) 
is also plotted against slot area ratio. It will be seen that 
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Fig. 8. 
Comparison of mass ratios with sonic and supersonic injection. 


(i) The ratio (supersonic injector mass flow)/(sonic injector mass flow) 
remains very nearly constant, to within about one per cent., over the 
range of « considered. 

(ii) Under the conditions assumed, supersonic injection leads to a reduction 
of about 13 per cent. in the rate of mass flow of injected air compared 
with that required for sonic injection. 


It is found that conclusion (i) holds for other values of M, and », for M, >1, 
whence it follows that a conyenient method of finding the mass ratios for sonic 
injection corresponding to finite values of o is as follows: — 
(a) For the given values of M, and », find the “ideal” mass ratio for super- 
sonic injection (M,=M,), Fig. 6. 

(b) Hence find the mass ratios corresponding to various values of « from 
equation (17). 

(c) Multiply these figures by the ratio (“ideal” mass ratio with sonic 
injection)/(“ ideal” mass ratio with supersonic injection), as given in 
Fig. 8, for the appropriate values of M, and 1. 

This method is considerably more convenient than the graphical method 
previously described. 


Fig. 6. 
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Fig. 8. 
Comparison of mass ratios with sonic and supersonic injection. 


(i) The ratio (supersonic injector mass flow)/(sonic injector mass flow): 
remains very nearly constant, to within about one per cent., over the 
range of « considered. 

(ii) Under the conditions assumed, supersonic injection leads to a reduction 
of about 13 per cent. in the rate of mass flow of injected air compared 
with that required for sonic injection. 


It is found that conclusion (i) holds for other values of M, and », for M, >1, 
whence it follows that a convenient method of finding the mass ratios for sonic 
injection corresponding to finite values of o is as follows: — 
(a) For the given values of M, and », find the “ideal” mass ratio for super- 
sonic injection (M,=M,), Fig. 6. 

(b) Hence find the mass ratios corresponding to various values of « from 
equation (17). 

(c) Multiply these figures by the ratio (“ideal” mass ratio with sonic 
injection)/(“ ideal” mass ratio with supersonic injection), as given in 
Fig. 8, for the appropriate values of M, and ». 

This method is considerably more convenient than the graphical method 
previously described. 
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Fig. 9. 
Comparison of the theoretical variation of mass ratio with injector slot area ratio with 
- values from actual wind tunnels. 


The curves of the relative injector mass flows under “ideal” conditions given 
in Fig. 8 show that supersonic injection leads to a decrease in the amount of injected 
air required, other things being equal, the saving being 20 per cent. or so compared 
with that required for sonic injection at reasonably high working section Mach 
numbers. While the values in Fig. 8 represent an increase in efficiency of the 
tunnel alone, the percentage increase in maximum running time for an induction 
tunnel driven from compressed air storage is found to be a little less than would 
be expected from a consideration of the “ideal” mass ratios. 
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6. Results and Conclusions 


Curves have been derived giving “ideal” mass ratios for induction-type high- 
speed wind tunnels over a convenient range of working section Mach number and 
diffuser efficiency. These curves are given for both sonic injection and supersonic 
injection (at the working section Mach number). A simple method is derived for 
obtaining mass ratios corresponding to finite injector slot area ratios from the ideal 
values, which correspond to negligibly small slot areas. From a consideration of 
these curves it has been shown that 


(i) For the tunnel alone, the highest mass ratio is obtained by using the 
smallest possible injector slot area and the greatest possible blowing 
pressure. 


(ii) For a tunnel driven from compressed air storage there is an optimum 
combination of slot area and blowing pressure, other things being 
constant, which gives the maximum running time”. 


(iii) There is no obvious reason why induction tunnels should not operate at 
appreciable higher Mach numbers than are attained in present practice, 
although the overall efficiency of the plant will decrease rapidly as the 
Mach number increases above about 2:0. Presumably the present diffi- 
culty in attaining high Mach numbers is due to factors such as non- 
uniformity of the stream at entry to the diffuser. 


(iv) The mass ratio shows some improvement when the injected air enters 
the mixing section at the working section Mach number. The increase in 
efficiency depends on the characteristics of a given plant, but should be 
of the order of 15 to 18 per cent. at reasonably high supersonic speeds. 


In practical induction-type high-speed wind tunnels, the mixing is not usually 
a constant-area process and the conditions of the stream are far from uniform at 
entry to the diffuser. A direct comparison between experiment and the theory given 
in this paper is therefore difficult. However, it will be seen from Ref. 4 that the 
general behaviour of a practical tunnel is similar to that predicted by the theory 
given. Fig. 9 shows that over a reasonable range of slot area ratios, the influence 
of slot area on mass ratio is closely in agreement with that predicted by the final part 
of the theory, although considerable discrepancies may occur when @ is less than 
about 0:04. However, it would seem reasonable to suppose that a practical tunnel 
would also behave approximately as predicted when the injection is supersonic. 


Under either of the injection conditions considered here, the static pressure of 
the injected stream at exit from the slot is greater than the local main stream static 
pressure, and hence further expansion of the injected stream occurs after it leaves 
the slot. Since mixing occurs over a considerable distance and, in fact, is never 
complete, initially the main stream suffers an effective contraction downstream of 
the slots, which may lead to “ choking.’ 


It is, of course, possible to obtain M,=1 
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as a solution to equation (8) by a suitable choice of the initial conditions of the two 
streams, but the theory given indicates that the normal running condition of a 
tunnel approaches this solution only when the working section Mach number is near 
unity and choking in the mixing section would be attained only by admitting a 
considerable excess of inducing air (“ over-blowing ”) at other working section Mach 
numbers. In practice, however, the effect of the expansion of the inducing stream 
after leaving the slot is considerable, and induction tunnels often operate near the 
choking condition when the working section Mach number is near 1:8 or so. To 
achieve even this value, it is usual to apply a considerable divergence to the tunnel 
walls in the injector, and it is not unlikely that this choking, due to insufficiently 
rapid mixing, provides the main limitation to the operation of induction tunnels at 
higher speeds. An analysis of this choking effect would involve, however, a detailed 
consideration of the mixing process. 
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A Note on the Performance of Ducted Fans 


BRYAN THWAITES 
(Imperial College) 


Summary: The strip-theory equations of flow through a ducted fan were the 
basis of a method"? of designing a fan for a given performance. The merit of the 
method was that it amounted to an explicit solution of the design problem, no 
iteration being necessary. 


In the inverse problem, of determining the flow past a given fan in conditions 
other than those of the design, two relations—the motor or tunnel characteristics, for 
example—will be given between the four quantities P, (2, , u (the motor power 
and angular velocity, the fan efficiency and the air speed at the fan). The problem 
is to determine the values of these four quantities together with details of the flow 
past any blade element. 


This paper shows that to a high degree of accuracy the performance of a ducted 
fan can be summed up by saying that P/{2* and » are both functions of {2/u only, 
functions which depend only on the fan shape. Solution by graphical means is 
therefore almost instantaneous. 


The fan-tunnel combination is shown to be governed by extremely simple 
equations which amount to saying that the tunnel power factor is a function of 
Q/u only which, again, can be calculated once for all for a given fan and tunnel. © 


Notation 
r radial distance y  tan~'(C)/C,) 
R, _ tip radius 6 blade angle 
dO n fan efficiency 
Q torque= | t Q/a 
duct 
R 
r root radius 
T  thrust= | dr R, 
r dP 
r, P  power=2Q= 
c chord of fan blade ae 
axial velocity © angular velocity 
a incidence N number of blades 
C, lift coefficient a, _lift-curve slope 


Paper received August 1952. 


{The Aeronautical Quarterly, Volume IV, February 1953) 


179 


| 


BRYAN THWAITES 


—a,  no-lift incidence 
Cy drag coefficient F tunnel power factor 
pair density expansion ratio of tunnel 


P,,P,,P, defined in equation (14): S,,5,,S.,S,,S, defined in equation (20). 
Q,’,Q,’,Q2.’ defined in equation (12): R,’,R,’,R,’ defined in equation (18) 


Note: A barred quantity equals the ratio of that quantity to its designed value. 


1. Strip Theory Equations 


The familiar form of the equations is given in equation (1) of Ref. 1 which 
may be re-arranged as 


= . . (i) 


tan g= 4 {14 cin + y)cosee sec sec } 


aT _ cot(+y)dQ 
dr r dr (4) 
tan y=Cy/C, . . (5) 


It is assumed that the fan is given in terms of @ and the section at any radial 
distance. 


2. The Expression for P 
Knowledge of the blade section implies the equation C,=a, sin(z+2,) which 
is, from (3), ; 
C, =a, sin (6 + a, —¢) ‘ ‘ (6) 


a, and @, are, of course, functions of r. Elimination of C,, first between (1) and (2), 
then between (2) and (6), gives 


8xr 
=sin (0 + 2% [ ] 
U=sin (6 + a,) tan y cot? Hoa. +sin (6+ 2, —y)sec y | cot cos (6+ 
(8) 


Elimination of cot ¢ would give a relation between dQ/dr, 9, u and the fan’s shape, 
but the variables could not be separated and no simple result could be obtained. 
Such a result is obtained at any rate for y=0, for then, from (8), ucot®@ is linear 
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in u and so, from (7), dQ/dr is a quadratic in u, the coefficients being known 
quantities. This suggests solving (8) for cot as an expansion in y. Let 


rQ 


— —cote= | (9) 
+ sin (0+ z,) 

Nea, 


(9) is the solution of (8) when y=0, and when y +0, substitution of this expression 
for cot ¢ in the terms of (8) other than those involving tan y gives 


K =cos (6+ z,) cot sin (6+ 2,) cot? (0) 


Now if y and (cot ¢—r{2/u) are both regarded as small quantities of the first order, 
the contribution of K to cot®# is of the second order. A first order approximation 
to K in (10) will therefore result in a third order error in the value of cot? and 
this is negligible. To this order of approximation, cot@ may be put equal to 
rQ/u in (10) to define K, which in turn defines cot¢ in (9), substitution of which 
in (7) gives 


dQ _ pNea, sin (6+ 2,)tan y+ sec y sin (6 + 2, —y)—u*r cos (6+ 2,) 


dr Nea, (11) 
1+ sin (6 + z,) 
If Ota - sin (6+ 2,)tan y_ 
1+ sin (0+ 2,) 
Q,’= r’ sin (6+ 2, —y)secy 
Nea 
1+ - “sin (8+ 2,) 
1+ sin (6+ z,) 
(11) can be rewritten 
R, 
Np 
and with | adr, n=0, 1,2, ‘ . (14) 


integration of (13) and multiplication by 2 give 
P=0°P, + uP, ‘ (5) 

This is the first principal result of the paper. 
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In the calculation of the Q,’, it is consistent with the approximation made to 
take constant value of y, typical of the blade section at a moderate C,; otherwise 
values of y may be taken from the design of the fan. The Q,’ should be tabulated 
during the computation, since they are required later for the calculation of C,. 


3. The Expression for 1 P 
From (4),* 
dT dPiu 


To the second order of the quantities y and (cot¢—r/u), with the value of cot ¢ 
from (7), ; 


(16) 


2-2 


dr 


? 
(17) becomes cot (+7) =——R,'+ —R,'+R,’ ‘ (19) 
r u u 
since dQ/dr in (17) is given by (13). 
R, 
With S,=4Np | ca,Q,’R,'dr 


R, 


S,=4Np | ca, (Q,’R,’—Q,'R,’) dr 


r, 


R 


S, = 4 Np Cay (Q.’R A = Q,’R.’) dr 


S,=4Np | ca,Q.’R,/dr 


*n here is the fan's efficiency not the (variable) blade element efficiency. Thus this equation 
should be interpreted as defining dP/dr whose integral from ra to R, is P. 


182 


rae Q,’Nca, 
With R,=rtan y+ 
‘Nea, 
R,'=1- 
p= 2: Nea _ tany 
r 
= 


PERFORMANCE OF DUCTED FANS 


we get from (16) and (19) 


This is the second principal result of the paper.* 


It should be noted that variation of u can be taken into account. If u=Uf(r), 
f(r) being a known function, the integrand of P,, has an additional factor [f (r)]", 
n=0, 1, 2, that of S, a factor of [f(r)]"-', n=0, 1, 2, 3, 4; while U replaces wu in 
(15) and (21). 


4. Method of Solution 


Equations (15) and (21) together with the two other relations between P, ©, », u, 
which must define the operating condition of the fan suffice for the calculation of 
these four quantities. Solution, especially by graphical means, is extremely simple 
since (15) and (21) amount to saying that each of P/9°, P/u’®, n, »P/u®, »P/Q® are 
functions of £2/u which can be calculated once for all for a given fan. The following 
are typical cases : — 


(i) P and {2 given (motor performance measured). P/{)* defines a value of 
©/u from which u and » follow. 

(ii) P given, and »P given as a function of u (motor output measured and 
duct power characteristics known). Since P/u* is a function of 2/u, u 
can be determined as a function of {2/u, and so therefore »)P/u’* also, 
which when compared with the known fan function for »P/u* gives a 
value for 2/u. u, © and » immediately follow. 

(iii) P given as a function of ©, and u given (motor control setting known and - 
duct velocity measured). The data can be put in terms of P/* as a 
function of 2/u, which when compared with the known fan function 
P/Q gives Q/u. ©, P and » all follow. 

(iv) u and »P given (duct performance measured). 1P/u* defines a value 
of 2/u from which ©, » and P follow. 


5. Blade Element Characteristics 


When u and 2) have been determined, dQ/dr may be found from (13). cot¢ 
and so ¢ is then found from (7), and so C; from (6). dT/dr is given by (4) and the 
pressure rise across the fan by (dT /dr)/(2tr). 


6. Fan in a Wind Tunnel 


This is the case to which the preceding theory will most often be applied, and 
particularly simple results are obtained. If the power factor of the tunnel is F, 


*It should perhaps be stated explicitly that equations (15) and (2 1) are solutions of ‘the strip- 
theory equations to O(7”), O [(cot ¢—r2/up] and O(dy/dC,). 
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TABLE I 
R=2-25 r,=1-189 N=3 
c=! @, =5-753 
r 6° Q,’ Q.’ 

1-189 | 1:0 0-01109 0-5158 0°8151 
1-366 21-9 1-1 0-01706 0-6235 1-0211 
1-543 19-0 1:2 0-02474 0:7275 1-2246 
1-719 16°8 0-03455 0-8300 1-4229 
1-896 14-9 1-4 0-04621 0-9236 1-6218 
2-073 13-4 0-06058 1:0165 1-8167 
2-250 12-1 1-6 


0-07615 1-1020 2:0104 


P,=0-0008199; P,=0-01790; P,=0-03095 
S, =0-00004205; S,=0-0001059 
S,=0-01752; S,=0:02429; $§,=0-005449 
Designed values: u=77 Q= 157 P= 8250 
n=0°892 F=0-222 e=2:303 


the fan duct area is A and equal to e times the working section area, then* 
nP=4 pu’ Ae*F 


Substitution of this value of P in (21) gives a quartic in 2/u, the value of which 
follows at once. This shows that the ratio 2./u for a fan-tunnel combination 
depends only on the power factor and the fan design, and not on the motor control 
setting. Substitution of this value of 2/u in (15) gives a constant value of P/0®, 
which again does not depend on the motor control setting. If this setting is given, 
the relation between P/* and © is therefore known from the motor characteristics 
and so 2 is immediately determined. u, P and 7 all follow at once. 


Solution of the fan-tunnel problem, given the motor control setting, is thus 
a matter of a few minutes. 


7. An Example 


The fan designed in Ref. 1 is defined by the first two columns of Table I. 
together with the values at the head of the table. Values of y are those of the design. 
The Q,’ are tabulated, for the calculation of the blade element details requires them. 
The R,,’ are not shown, but the final coefficients P, and S, are given at the foot 
of the table. 


*Contrary to usual practice, the fan output is taken in this formula, for then there is more 
justification in regarding F as a constant, independent of u; but there would be no difficulty 
in adopting the usual definition. 
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It is convenient to work in terms of the ratio of a quantity to its design value*, 
which will be denoted by a bar. Thus, P=8250 P. If t=Q/u, equations (15) 
and (21) may be written 


4:0759 3-4566 
t 


p-0 (0-3807 + (22) 


F= = —0-04476 t* — 0-05529 + 4-48591 — 3-05026 t—0-33560. (23) 


These equations are, of course, satisfied if all quantities take the value unity. 
Graphs of P/Q’, 7P/u*® and F may be drawn against t, and solution is a matter 
of minutes. 
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APPENDIX 


In Ref. 1, the performance problem was dealt with briefly by an exact solution 
for ¢ from equation (8), which of course prevented the simple results of the present 
paper being obtained. However, the solution for ¢ given by equations (18) and (19) 
of the original paper is incorrect, and these equations should therefore be 
replaced by 


rQ 
u Nea, 


Nea, 


cos y + cos (6+ 2, +7) 


_ 


cos y+ sin (0+ 2,—y) 


rQ 
u Nea, 


in (6 


cos y + cos (8+ a, 


tan d= 


The form for ¢ given in the present paper’s equations (9) and (10) is however 
far preferable. 


*Or to any fixed value for which the solution has been worked out. 

¢In fact, the right-hand sides of the two equations have been multipled by factors 0:9898 and 
0-9927 respectively. The differences between these factors and unity represent the error 
arising out of the theory. 
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On the Resistance of Screens 


K. E. G. WIEGHARDT 


(Formerly Admiralty Research Laboratory, Teddington, 
now at Hamburg University) 


Summary: The resistance of various screens in a tube is reduced to the drag of 
a single piece of wire by assuming that the average velocity in the screen itself rather 
than the free stream velocity is the significant factor. Almost all the available tests 
on screens or gauzes are brought into line with a single experimental function similar 
to the drag of a cylinder plotted against Reynolds number. 


Notation 
wire diameter 


width of the square mesh 

porosity =(area of holes)/(total area) 
velocity before and after the screen 

static pressure before and after the screen 
density of fluid 

kinematic viscosity of fluid 

Pi-P2 


p 


resistant coefficient = 


B? 
i-p* 


Cp drag coefficient of an infinite cylinder 


c reduced resistance coefficient = 


1. Normal Flow Through Screens 

The resistance of screens or gauzes in a tube has been measured several times 
since 1932°-, Also some theories were established, two of which':*) lead to 
formulae deriving the resistance only from the geometrical dimensions. Considering 
the complexity of the general problem of aerodynamic resistance, they are likely to 
be only rule-of-thumb formulae, even if they cover the tests in a certain range, and 
so a less ambitious approach might be preferable. 


Attention is focused only on screens made of round wires forming a square 
mesh, as these are defined geometrically alone by the wire diameter d and the 


porosity 


area of notes _ ( 


] ‘ ‘ (1) 


area 
where / denotes the width of the square mesh. 


Originally received March 1952; revised August 1952. 
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The resistance coefficient k is defined as the pressure drop through the screen 
divided by the pitot head in the straight tube, i.e. 


(2) 


The proposal has been made"? to reduce the resistance of the screen to that 
of a single wire, i.e. an infinite cylinder. As a fraction of (1-8) of the whole area 
is covered by pieces of wire with a drag coefficient c, the resistance of the screen 
is expected to be roughly equal to the sum of the drag of all these wire pieces, or 


k=(1-B)c. . ‘ ‘ ‘ (3) 


In the same way as c depends on the Reynolds number u,d/v, k/(1 - 8) should 
be a function of u,d/v alone for all screens. 


The experimental results, since they give various curves for various values of 
8, fail to substantiate this. The reason appears to be that the characteristic velocity 
for the flow round the wires is not the one in front of the screen u,, but the higher 
average velocity through the screen, u,/8. This idea has also been used by 
G. I. Taylor™. 


Hence, for computing the drag it may be assumed that the single wires are in 
a field of flow with free stream velocity u,/8. This may be a bold assumption 
because the actual velocity in front of the wires is u,; on the other hand, it is the 
velocity (outside the boundary layer) on the rear part of the cylinder which 
influences the wake and therefore the drag. With this assumption the resistance per 
unit area becomes 


=K 


Now c should depend only on u,d/(fv). That this is so is shown in Figs. 1 and 2, 
where all the available tests are plotted on logarithmic scales. Although k varies 
between 0-35 and 19, all the measurements lie along a single curve, particularly the 
very accurate tests from Ref. 6 where, incidentally, an empirical formulae 


is given for high Reynolds numbers. 


Of this whole series from Ref. 6 only the points marked (O) for a single screen 
are out of line, but in this case the mesh was rectangular with sides in the ratio 1:1-25, 
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The reduced resistance of gauzes, from Refs. 1, 3, 5, 6. 


while in all other cases the mesh was square; the discrepancy presumably indicates 
the influence of the geometry of the system. Other points which do not fit well are 
the highest values of Ref. 1 (marked x). No explanation can be offered here, and 
it might be useful to make some new tests at high Reynolds numbers, e.g. in a water 
tunnel. In a water tunnel with a large working section it would be possible to 
achieve even turbulent flow round the wires (critical Reynolds number u,d/(8¥) 
about 3 x 10°). 


Further, it is seen from Fig. 1 that the difference between the average c and the 
drag coefficient cp of a single cylinder in infinite flow is unexpectedly small, in view 
of the crudeness of the whole reduction. In particular no attempt has been made to 
take account of the fact that the wires composing a mesh are not straight, nor of the 
influence of the points at which wires cross one another. 


2. Comparison with Previous Theories 


Another theoretical relationship for k has been derived in Ref. 5. Apparently 
the main interest lay in the resistance of perforated sheets. Hence the flow through 
such a sheet or a screen was reduced to that through a single hole as at the vena 
contracta. This gives 

k= BC —1, ‘ (7) 
with C=contraction coefficient and y=fraction of the pressure loss which is regained 
when the stream again becomes uniform behind the sheet. Theoretically this 
approach is as good as the previous one; but when considerations are restricted to 
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Fig. 2. 
The reduced resistance of gauzes, from Ref. 7. 


screens made of round wire it seems more practical to take as basis the drag of a 
cylinder than to find out how C and y depend on 8 and the Reynolds number. 
Furthermore, with most screens, when £ is not too small, the mesh length or the 
distance between the wires is several times the wire diameter. Therefore the com- 
parison with the flow around a cylinder in an infinite fluid seems more natural than 
that with the jet through a single hole. For screens, it was found that for high 
Reynolds numbers 


0-70 


Assuming an average value of u,d/v of 200, this relation gives the curve “T.D : Th.” 
in Fig. 1. (For other values of u,d/v the same curve has to be moved horizontally.) 


From a paper of Eckert and Pfliiger a graph of k against 8 for u,d/v=200 
is also plotted in Fig. 1, giving the curves “E.P : Exp.” and “E.P: Th.” Their 
theoretical formulae is ; 

_ (1-8); 
k= (45°) 


In Fig. 1 this curve “E.P : Th.” does not fit the tests at all whereas in their 
plotting the agreement is reasonable, though not very good; a striking example of 
how much depends on the art of presentation. 


Glauert, Hirst and Hartshorn” derived the formula 


3). 


which gives the curve “G.H.H: Th.” in Fig. 1 again assuming u,d/v=200. 
189 


LASS 


K. E. G. WIEGHARDT 


y 


| 


10 
100 
100. 


ie 
0. 


4 


Fig. 3. 
Nomogram for the pressure drop ina gauze in air flow. 


Figure 1 suggests another formula. For u,d/(8v) between 60 and 600 all the 
tests at various Reynolds numbers and porosities can be approximated by the broken 


line 


On the other hand, the American tests in Fig. 2 give lower values, 


No explanation for this difference is obvious. Admittedly the test arrangements 
were not quite the same. At the National Physical Laboratory the screens were put 
into a tube, while at the National Bureau of Standards the screens were fixed at the 
end of a square duct discharging into free air. Still, the flow through the screen 
should have been the same. As there might have been a difference in the surface 
roughness of the wires forming the gauzes, it would have been useful to exchange 
the gauzes of the two test stations. Incidentally, the dotted line K in Fig. 2 
corresponds to a silk cloth manufactured for bolting flour. 
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Summarising, an average of all the tests may be assumed from Figs. 1 and 2 
to be - 


or 


This formula has been used in the simple nomogram in Fig. 3 for calculating the 
pressure drop of a gauze in air flow. 


For example, if d=; in., 8=4, u,=20 ft./sec. 


ud ~ 350, 
c= 0-78, 


and the pressure drop in the gauze is 


| 
Pi—Pa= 5," x xC 


=0-1 x2 x0°:78=0-16 in. water. 


3. Critical Reynolds Number for Eddy Production 


Schubauer, Spangenberg and Klebanoff'” also observed by a hot wire the eddy 
shedding immediately after the screen and gave a graph for the critical Reynolds 
number, for the beginning of eddy shedding or of the initial production of turbulence, 
against 8. Obviously a straight line is a good approximation. Yet this gives a 
-constant critical Reynolds number if it is referred again to the velocity u,/f in 
the screen: 


(42) 


A Reynolds number of the same order is the smallest at which a vortex street 
is formed past a single cylinder in infinite fluid as is to be seen from F. Homann’s 
photographs of the wake“. 
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Flow Changes in Gases in which Mass 
and Impulse are Conserved 


L. G. DAWSON 
(Rolls-Royce, Derby) 


Summary: This note discusses changes in the state of a gas flowing in a duct 
of constant area. In the past, changes, such as normal shock waves, combustion 
and other phenomena, which are defined by the equations of energy, conservation 
of momentum and mass flow, have each been treated on their merits. In this note 
a method is developed whereby all phenomena governed by these three equations 
can be solved by a single general method. The method rests on the derivation of a 
parameter which is unaltered in value by the change, in all cases where the total 
temperature is constant. A shock wave is an example of such a discontinuity. In 
problems of heat addition or extraction, the parameter changes its value only because 
of the change in total temperature. The change in total temperature may be 
calculated from the known quantity of heat added or extracted. 


The parameters derived are useful in showing how problems of this type should 


be attacked analytically. With oblique waves it is easy to derive a relation between 
the normal velocities before and after the wave, and it is probable that this 
relationship has not been published before. 


1. Introduction 


There are many types of change in gas flow problems which are governed by 
the same laws. The change from supersonic to. subsonic flow through a normal 
shock wave is described by the same physical equations as combustion in a 
parallel tube in a subsonic stream. In the past these different processes have each 
been treated on their merits, but it should be possible to treat all flows in which the 
changes are governed by the same equations in a perfectly general manner. The 
purpose of this note is to show that all flow changes which are specified by 


(i) conservation of mass flow, 
(ii) an equation involving change of momentum, 
(iii) the equation relating the total energies before and after the change 
may be described in a single general treatment. 
The usual solution of a particular problem is to write down the three equations 


and to solve for the particular variables of interest. Graphical methods have been 
evolved for solving the equations for a particular type of process; for example, 
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curves may be drawn showing the velocities before and after a normal shock wave, 
but no general method has been produced for treating all the eed which are 
governed by the three equations already mentioned. 


The object of this note is to demonstrate a general method which will supply 
all the information normally obtained by solving the three equations for the 
variables required. A great advantage of this method is that it provides a mental 
picture of the results of the three equations taken in conjunction and often shows 
how a particular calculation should be carried out in principle. 


If, instead of solving the three equations, a parameter is derived which contains 
numbers which represent the mass flow, the change in momentum and the total 
temperature, then this parameter will remain constant in any change in which the 
total temperature is constant. When there is a change in total temperature, this 
will be known; for example, in a combustion problem, it will be defined by the 
heat released. 


The new approach will be clearer after the parameter has been derived. 


Throughout the calculations, the ratio of the specific heats has been considered 
to be constant and to have a value of 1:4. In some of the problems considered 
there is a change in the ratio of the specific heats, particularly in combustion 
problems. The curves which have been calculated are intended to help to convey 


the ideas involved; the correct value of y can be used for particular cases. 


Notation 
Suffixes “ 1” and “2” denote conditions before and after a change. 
a__ local speed of sound 
speed of sound corresponding to stagnation temperature 
area normal to flow 
impulse= A (p + pu’) 
static pressure 
total pressure 
static temperature 
total temperature 
velocity normal to wave 
velocity parallel to wave 
absolute velocity 
mass flow 


cw Hn FS 


stream density 


specific heat at constant pressure 
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ratio of specific heats 

gas constant 

wave angle 

angle through which flow is deflected (Fig. 5) 


Derivation of the New Parameter 


The mass flow is defined by the equation 


W=pAu. 
Using the energy equation 
K,t+4u?=K,T 


and the speed of sound corresponding to stagnation temperature 


these two equations give 


and the mass flow equation may be written 


The impulse of a flow is defined as 


and this quantity is.conserved if there are no external forces acting on the system. 


This may be written 


Y 
| a 
B 
a,?=yRT, 
¥ 2 a?’ 
| | | 
yw T 
Ap 
2 
I 
giving Ap 
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and it is more convenient to deal with the reciprocal 


i 
y+1 uv? 


1+ 


Multiplying equation (1) by equation (2) gives 


|(R), 
Ap Y I yt 


W/T (2) 1 
I y/ 
1+ a 


The left hand side of equation (3) is the parameter required. 


3. Use of the Parameter 


As an example, consider conditions before and after a normal shock wave. 
The mass flow per unit area is the same, and the change of momentum in the 
direction of flow is solely due to the gas pressures. Thus W is constant, and 
I= p-+pu? is constant before and after the wave. The stagnation temperature is 
unchanged, making the left hand side of equation (3) constant. 


If it is the changes in total pressure which are of interest, equation (2) is written 


which is the most useful form. Again the left hand side is plotted and changes 
only by virtue of the changes in total pressure. 


In Fig. 1, (W/T/I)/(R/y) has been plotted as a function of Ap/I and 1/(AP), 
the points being labelled with the appropriate values of u/a,. 


In Fig. 2 (W/T/I)/(R/y) is plotted as a function of u/a,. 
196 


or 
(3) 
a3 P 
a3 
1 = 
But 
L 2 a,” 
1 


FLOW CHANGES IN GASES 


\ 


al 
* 
| 


| 


o6 o7 ose o9 
Fig. 1. 
Impulse and continuity function for y= 1-4. 


At this stage in the discussion, only one of the curves is significant, that is, 
the curve for which the parameter v/a, is zero. The significance of the other curves 
in Figs. 2 and 3 will be clear after the discussion of the oblique wave case. 


The use of the parameters which have been defined so far is best made clear 
with four examples. The numerical values of the parameters for each example are 
shown in Table I, and a short discussion of each type of change follows. 


3.1. Normal shock 


The approach velocity u/a, has been chosen as 1-4, which is supersonic. Since 
it is a normal shock there is no change in total temperature. (W/T/D/(R/y) 
may be found from either Fig. 1 or Fig. 2 and is 0-418. This number is unchanged 
through a shock wave and this enables the remaining quantities to be determined 
from Fig. 1. The values of Ap/J at approach and on leaving are determined, the 
ratio giving the static pressure ratio across the system; this is 3-94, the static pressure 
behind the wave being higher than ahead. 

The loss in total pressure is determined from the values of //(AP), the ratio 
of this parameter giving the ratio of total pressures behind and before. This ratio 
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4 


do 


Fig. 2. 
Continuity function. Oblique waves. 


is always less than 1:0. The leaving velocity is now subsonic at a value of 
u/a, of 0-595. 


3.2. Normal shock with heat addition 


This system has not occurred experimentally with the wave stationary relative 
to the observer, but it is the stationary wave counterpart of a detonation wave. The 
approach velocity has been chosen to be the same as in the first example and the 
total temperature ratio has been chosen to be the maximum possible. Inspection 
of Fig. 1 will show that the maximum value of (W/T/D/(R/y) is 0-456, so that 
the total temperature ratio is (0-456/0-418) = 1-19. 


The static and total pressure ratios are determined from the appropriate values 
of Ap/I and I/(AP) as in the first example. The velocity leaving the system is 
now sonic because the maximum amount of heat has been added. 
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Fig. 3. 
Impulse function. Oblique waves. 


3.3. Subsonic combustion 


This is a far more practical example. A low value of the appropriate velocity 
has been chosen and an arbitrary total temperature ratio of 4-0. Again the value 
of (W./T/I)/(R/y) changes only by virtue of the change in total temperature. 
The value of the parameter corresponding to the approach velocity ratio u/a,=0°1 
was 0-095, so that after combustion the parameter has the value of 0-095 x 2=0-19. 
This allows the leaving velocity ratio to be found and defines the values of Ap/I 
and [/(AP) which determine the pressure ratio. 


3.4. Limiting combustion temperature in a parallel tube 


The amount of heat that can be added in subsonic combustion is limited and 
can be determined by the temperature ratio required to give a maximum possible 
value of (W/T/I)/(R/y) which corresponds with sonic leaving velocity. Here 
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the temperature ratio is not known, but is determined from the limiting value of 
(W/T/D/(R/y), which is 0-456, and the value for combustion, which in this 
example is 0:275. The limiting total temperature ratio is then the square of the 
ratio of these values, i.e. 2°75. 


4. Oblique Waves 


Once the principle suggested here has been grasped, it can be extended in a 
variety of ways. It may be used to define the changes which occur in an oblique 
wave. The parameter (W/T/I)\/(R/y) is generalised so that the flow W is defined 
only by the velocity normal to the wave u/a, and is independent of the velocity 
parallel to the wave v/a,. The total temperature has the usual significance but / 
is defined as p+pu?, where p is the static pressure of the stream and p the static 
density; u is the velocity component normal to the wave. When the detailed 
expression for (W/T/1)/(R/y) is calculated it will contain both the velocity normal 
and parallel to the oblique wave. The derivation is quite straightforward and will 
be found in the Appendix, together with the expression for //(AP). In Figs. 2 
and 3 the general expressions for these two parameters have been plotted with 
v/a, as a parameter. 


The way in which the continuity function is used for oblique waves is not quite 
the same as in the plain wave problem. If a constant value of the continuity 
function (W./T/I)/(R/y) is chosen, then the value of the velocities normal and 
parallel to the wave are related and the one may be plotted in terms of the other. 
This is shown in Fig. 4 where the horizontal axis is u/a, and the vertical axis is v/a). 
Any one curve defines the relation between these velocities for a fixed value of the 
continuity function. 


It will be noticed that for a fixed value of the continuity function and v/a, 
there are always two possible values of u/a,. This is simply a description of the 
fact that it is only the velocity normal to an oblique wave which changes. The 
line joining the minima of the curves corresponds to the velocity normal to the 
wave being sonic. There will be no discontinuity if the normal velocity is less than 
sonic. If the wave angle of an oblique wave is z as defined in Fig. 5, then a line 
drawn in Fig. 4 making an angle 2 with the vertical axis represents the absolute 
velocity, provided that it is drawn to the correct length. The end point of this 
line will always lie on a curve corresponding to some constant value of the continuity 
function. Only a small number of these curves have actually been drawn in Fig. 4 
but they can be interpolated for any particular case. 


So long as the magnitude and angle of the absolute velocity correspond to a 
normal velocity greater than sonic, then there will be a smaller alternative value of 
u/a, which is defined by the point where the curve corresponding to the constant 
value of the continuity function shows the same value of v/a,, since the velocity 
parallel to the wave cannot alter. 
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As an example, consider a wave at 45° to the axis, with absolute velocity / 2a,. 
This intersects the curve of the continuity function of value 0-5 at u/a, and v/a, 
equal to 1:0. The outlet condition from this wave is shown, and corresponds to 
u/a,=0-66, the stream being deflected through an angle of 114°. The diagram may 
not have great virtues for solving this simple problem, particularly as it will be 
shown that by considering the behaviour of the continuity parameter, it is possible to 
derive a simple analytical expression for the velocities before and after the oblique 
wave. However, it gives a solution of a more complex discontinuity. If heat is being 
extracted along an oblique line, for example, by a condensation process, then the 
total temperature is no longer constant and the continuity function must be reduced 
in proportion to the square root of the change in total temperature. If the total 
temperature ratio is 0-9, then in the example already considered, the continuity 
function would change from 0-5 to 0-475, and the outlet condition would be defined 
by u/a,=0-555, the deflection now being increased to 16°. It will be interesting to 
know if this phenomenon can occur physically. 


/ 


/ 


Fig. 4. 
Oblique wave diagram, 
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Fig. 5. 
(a) Oblique wave. (b) Oblique wave diagram skewed 
to fit velocities. 


By equating the values of the impulse function for the oblique wave, a simple 
relationship between the normal velocities before and after the wave, is obtained : — 


aa, ytil 


This is a generalisation of the well-known relationship for normal waves :— 


With the aid of the concept of the continuity function, the derivation of the 
more general expression is extremely simple, as shown in the Appendix. It is not 
known whether this has been done previously. 


APPENDIX 


The continuity function for an oblique wave is obtained from 


and therefore = 
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The impulse normal to the wave 


Since V?=u? +’, 


Dividing [W /T/(Ap)]/(R/y) by the first alternative gives 


The momentum equation parallel to the wave shows that v/a, is constant. 


To obtain an analytical relation between the velocities before and after the 
wave, the alternative expressions for (W/T/I)/(R/y) are equated : — 


Multiplying out gives 


sil pA t yRT | 
<i 
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The Arithmetic of Field Equations 


A. THOM, M.A., Ph.D., D.Sc. 


(Professor of Engineering Science, Oxford University) 


Summary: The paper describes in detail an older method than Relaxation of 
approximating to the solution of equations of the Laplace and Poisson type. The 
corresponding fourth order equations are discussed briefly, and a method of dealing with 
certain non-linear equations is indicated. A description is also given of the propagation 
of errors in the fields due to various causes. 


1. Introduction 


In the treatment of the subject given here it is assumed that the reader has no 
previous acquaintance with the method, but has a general knowledge of the differ- 
entiation and integration of functions given by a numerical table with equal intervals 
of the argument. Rigid proofs are not given, as these would make the paper much 
too long, but the method of proof is usually indicated. 


It should be made clear at the outset that the process is different from 
Southwell’s Relaxation Method. Since 1928 it has been applied successfully to 
hundreds of problems, some of which were of considerable complexity. The method 


is described by means of examples in the early part of the paper. The use of large 
squares or groups of squares is then explained and applied to various problems. 
This was published by the author in 1933, not knowing that, as far as the application 
to Laplace’s equation is concerned, it had already been used by others. 


Recently the method of groups of squares has again been mentioned by one or 
two writers, but the application to the double-field system of solving equations of 
the viscous flow type, advocated by the author in 1933")*’, had apparently been 
entirely overlooked until Woods drew attention to its advantages’. 


The indirect method of operating in the ¢,¥-plane was outlined in Refs. 1, 4 
and is fully illustrated here. The reader who wishes to look into the solution of 
non-linear equations will find that the methods given here are applicable, but 
reference should be made to various published papers for more complete details. 
' (See list of references.) 


Notation 


2m diagonal of square 
2n side of square 
r,9 polar co-ordinates in Section 3 
Received November 1952. 
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=x+iy 
see Section 5:2 


the value at any point in a square field due to A=1 at the centre 
square 


polar co-ordinates in w-field 

defined in Section 5 

stream and potential functions for grid (i.e. Laplacean Solution) 
velocity vector 

vorticity 

distance from wall in x, y-plane 

distance from wall in 7-plane 

maximum slope of side of swelling in wall 

maximum height of swelling 


in Section 8.1, the dash is used when needed to distinguish the 
disturbed values 


defined by Y=Y+p’=¥+yp, yw 
(log g)/y 
standard deviation 


General 


The problem is to find numerical solutions to equations of the type 


but attention will first be focused on the case V*~=0. The boundary values may 
be known but this is not always so and in two of the examples worked these have 
to be found as the solution proceeds. 


The method in its basic form depends on obtaining the value of ¥ at the centre 
of a square or diamond, called hereafter the unit square, in terms of the corner values 
and the differential equation. Take the diamond C,C,C,C, with diagonal 2m and 
side 2n (Fig. 1(a)). Express by Taylor’s theorem the value of ¥ at each corner (¥-) 
in terms of ¥ the value at the centre. Adding the four values gives 


in m*. . ‘ (2) 


If the diamond is not too large, i.e. if m is small, the higher derivatives can be 
neglected and 


4m? VY 


where ¥y is the mean of the corner values. Note that this is independent of the 
orientation of the diamond and so applies to any square. 
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Fig. 1(a). Fig. 1(b). 


3. The Case V*¥=0 


This leads to the very simple form 


If approximately correct values (or for that matter, any values) are written at the 
corners of overlapping squares or diamonds throughout a field where the boundary 
values are known and equation (4) is applied systematically to point after point, by 
sufficient repetition values will be obtained getting nearer and nearer to the solution. 
This is a laborious process but anyone who tries it seriously soon finds short cuts. 
If each square is to be treated individually the following routine seems to be the 
quickest method of approaching the settled field. It is recommended that the first 
two examples be actually worked step by step, otherwise the explanation is difficult - 
to follow. 


In Fig. 2 the values round the edges of the large square are to be assumed as 
fixed boundary values. They were calculated from 


Y= v/reos 5 ; (5) 


where r,@ refer to an origin just outside the lower left corner. This expression 
satisfies V?4=0 and so if (4) be used to fill in the field a check can be obtained 
from (5). 


A start is made by writing guessed values at all grid points, i.e. at the corners 
of all squares. The values used will have no effect on the final result (but see 
Section 9); nevertheless the better the guesses the sooner the solution is obtained"°’. 
To illustrate, suppose the values marked G have been guessed. These are written 
in. Assume now that it is decided to operate on diamonds of the size shown. 
Starting at A, write under the guessed value the mean, 2°14, of the four surrounding 
values. Do the same at B but in forming the mean use 2°14 for the value at A, 
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The errors which may be present due to the neglected terms in (4) must also be 
considered. Reducing the size of the working square will, it is true, reduce this 
error. Thus in Fig. 2 squares could be used instead of diamonds by inserting in the 
centre of each square on the sheet the mean of the corner values, and then working 
systematically over the whole sheet again. The process is now more than twice as 
long, since there are now 25 internal points instead of 9. Fortunately there are 
alternative procedures, one of which will now be described. 


4. Groups of Squares 


The whole field may be worked in a different manner by using groups of 
squares. This is more accurate and faster than subdividing, for reasons to be 
explained later. 


Consider the group of four squares shown (Fig. 1(b)). Assuming that the values 
C,A,C,A, and so on are known there are five unknown internal values, namely the 
central value ¥ and the value at each cross. An application of (4) to each of these 
points gives five equations which solve easily“ to give the value of wo in terms of 
the outside values as 


One advantage of this type of formula is that the actual values at the crosses need 
never be found or written in. The advantage in accuracy of a subdivision is gained, 
without actually subdividing. It is also a neater solution than that of filling in every 
square and is not so much affected by the end figure rounding. 


A different approach to the same idea is to express by Taylor’s theorem the 
values at the eight boundary points (Fig. 1(b)) in terms of the central values. 
Bickley’ shows that this leads to the following expression for a Laplacean field 


This gives a weight of 4 to the A values instead of the 2 given by (6). In Ref. 13, 
Woods shows the theoretical superiority of (7) in the general case, but particular 
cases may arise where (6) will give a better answer. Designating for convenience 
these formulae as the “twelve” and the “twenty,” the following rule has been 
found to be a reasonably good guide. It is based on a considerable experience of 
the two formulae. 


When a singularity or other similar disturbance lies near the group, use the 
“twelve” formula if the singularity is nearer an A than a C, and the “twenty” if 
itis neareraC. The relative weights given to the As and Cs by the two formulae 
give some basis for the rule. It is however inadvisable to change formulae on one 
sheet, and in the next example the “twenty” will be used throughout since it is 
slightly easier in practice. 
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Numbers in brackets ( ) show the order in which the figures were inserted. 


Repeat for C and proceed to the next line and the next again, always using the most 
up-to-date values in forming the mean. Returning to A, the mean using the new 
values might again be written, but there is no need to do this. Instead write at A the 
mean (which is 5 in units of the last place) of the differences or increments which 
have occurred at the surrounding points. Then at B write the mean of this 5 and the 
differences at C and E. 


Thus the sheet is covered systematically. At each point it is found that four 
differences are being used, all of which are new since the last operation here. Thus 
the solution is proceeding as rapidly as possible. If a difference works out as zero, 
record a zero and pass on; on the next round it is possible that a finite difference will 
appear due to a disturbance propagating itself from elsewhere. When ultimately 
all points carry recorded zeros, at each point the sum of the increments or differences 
is added to the first calculated value and so the final value is obtained. It is advis- 
able to repeat the whole process, using the values just obtained. This acts as a 
check but also reduces errors which may have arisen from rounding off the 
end figures, 
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The errors which may be present due to the neglected terms in (4) must also be 
considered. Reducing the size of the working square will, it is true, reduce this 
error. Thus in Fig. 2 squares could be used instead of diamonds by inserting in the 
centre of each square on the sheet the mean of the corner values, and then working 
systematically over the whole sheet again. The process is now more than twice as 
long, since there are now 25 internal points instead of 9. Fortunately there are 
alternative procedures, one of which will now be described. 


4. Groups of Squares 


The whole field may be worked in a different manner by using groups of 
squares. This is more accurate and faster than subdividing, for reasons to be 
explained later. 


Consider the group of four squares shown (Fig. 1(b)). Assuming that the values 
C,A,C,A, and so on are known there are five unknown internal values, namely the 
central value Wp and the value at each cross. An application of (4) to each of these 
points gives five equations which solve easily) to give the value of ¥. in terms of 
the outside values as 


One advantage of this type of formula is that the actual values at the crosses need 
never be found or written in. The advantage in accuracy of a subdivision is gained, 
without actually subdividing. It is also a neater solution than that of filling in every 
square and is not so much affected by the end figure rounding. 


A different approach to the same idea is to express by Taylor’s theorem the 
values at the eight boundary points (Fig. 1(b)) in terms of the central values. 
Bickley’ shows that this leads to the following expression for a Laplacean field 


This gives a weight of 4 to the A values instead of the 2 given by (6). In Ref. 13, 
Woods shows the theoretical superiority of (7) in the general case, but particular 
cases may arise where (6) will give a better answer. Designating for convenience 
these formulae as the “twelve” and the “twenty,” the following rule has been 
found to be a reasonably good guide. It is based on a considerable experience of 
the two formulae. 


When a singularity or other similar disturbance lies near the group, use the 
“twelve” formula if the singularity is nearer an A than a C, and the “twenty” if 
itis nearera C. The relative weights given to the As and Cs by the two formulae 
give some basis for the rule. It is however inadvisable to change formulae on one 
sheet, and in the next example the “twenty” will be used throughout since it is 
slightly easier in practice. 
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4.1. Example of the use of the “ twenty” formula 


Take the example already worked in Fig. 2 but increase by two the number 
of significant figures. The boundary values then become as shown in Fig. 3. As 
starting values for the 9 inside points the final figures from Fig. 2 are inserted. 
Since the adjacent “difficult” point, the origin, is near the corner of the adjacent 
group, the “twenty” formula is probably the best. Accordingly this is applied to 
each point in the same order as was previously used. The adjustments are 
made in the same order, using the “twenty” formula. It should be noted 
that although it appears that 9 internal points are being used, rather better accuracy 
is being obtained than would be given by the 25 necessary with equation (4). It is 
of interest to compare the result at the centre, 2:3869, with the true value from (5), 
namely 2:3878. Using the “twelve” formula in the same way gives 2-3853, 
showing as expected for this case a slight gain in accuracy by using the “twenty.” 


4.2. The “hundred” formula 


The effect of grouping four squares together has been studied. In exactly 
the same manner sixteen squares can be grouped (Fig. 4). Building up from the 
unit square, as in deducing the “twelve” formula, it is found? that 


476 ¥o=9 ‘ ‘ (8) 


Applied to the previous problem this gives the same result, in one operation, as 
that given by (4) applied to 25 points, or (6) applied to 9 points, but it gives the 
result without introducing end figure errors. 


A corresponding formula can be deduced for this group by the method used 
by Bickley in deriving the “twenty” formula, i.e. by using the 16 equations given 
by the application of Taylor’s theorem to the 16 boundary points. The symmetry 
of the expression allows the elimination of all derivatives and cross-derivatives up 
to and including the eleventh order and gives 


7548 vo = 83 +512 +780 ZA ‘ (9) 


It ought to be mentioned that in deriving (9) Laplace’s equation, V*/=0, has been 
used, and also in effect all the required derivatives such as 


_ 
* 


Thus in the form given (9) is limited to Laplace fields (but see Section 8). 


(10) 


The remarks which were made about the relative accuracy of the “twelve” 
and “twenty ” formulae in Section 4 apply equally to the comparison of (8) and (9). 


An application of (9) to the problem in Fig. 3 gives 2:3870 for the central 
value, while (8) gives 2:3853. The similarity of the results, in spite of the widely 


211 


6 
9 | 
8B 


Formula 


wert sin %,9 


Fig. 6. 
Corrections to be applied when using the “ twelve * or “ twenty ” formulae near a corner. 


different weights applied to the As, Bs and Cs by these two formulae, suggested that 
there might be a form simpler than either which would not differ much from (9). 
It appears that the following formula (the “hundred ”) satisfies these conditions. 


The fact that 100=4+7x8+10x4 makes this a likely form and on examination 
it proves to be almost as good as (9). In the example it also gives 2:3870. It will 
be found that the agreement with (9) is usually so close that the latter can be 
dispensed with and the “hundred” used instead. It is difficult to imagine a 
simpler formula, or one which, while bettering the accuracy of a sixteen squares 
subdivision, is so rapid in use. 


The “hundred” formula is particularly useful in working long fields four 
squares wide. Very little labour is involved in filling in the values on the centre 
line by the “hundred” and on the other two lines by the “twelve” or “twenty.” 
An example is given later. 
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In cases where these group formulae are not being used on the actual working 
of a field they form a very useful check. It is important to remember that all of 
them can be used in the diamond form, i.e. with the diagonal vertical. For example, 
in Fig. 5 the “twenty ” form might be applied to the diamond shown. To get to the 
line of points (along x, x) adjacent to a boundary (BB) without filling in additional 
internal points, a modified “twelve” is used“: — 


This necessitates extra values on the boundary. 


It might be mentioned that formulae intended to give values at points other than 
the centre of these large squares are not a success. Their essential lack of symmetry 
makes them awkward to use. It is much simpler to fill in the other points by the 
“twenty” formula. 


5. The Use of the Increment 4 
In many problems V?y is not zero. Then for the unit square 


(13) 


where A= 


In general, 4 varies from point to point and in non-linear solutions has to be found 
from the values or their derivatives as the solution proceeds. In working over a 
field with (13) the A is applied only once. It is of course not used in the subsequent 
application of differences, but in the case of non-linear equations it must be recalcu- - 
lated before the field is gone over again. Equations (6) and (7) can also be modified 
by applying a A. 


5.1. The reciprocal theorem 


Suppose that in a field squares of definite size are being used to operate on a 
function ¢ with fixed boundary values. Consider any two of the squares, one 
centred on S, and one on S,. Then the effect at S, of a A of amount A, applied to 
the square at S, is the same as the effect at §, of 4, applied to the square at S,. If 
definite groups of squares (e.g. the “twenty” formula) are being used the theorem 
also applies. Use will be made of this later. 


5.2. The application to Laplace’s equation 


Consider how the idea of using the increment A can help in the solution of 
V*~=0. Difficulty is always experienced when a singularity occurs on the edge of 
a field, or just outside the edge. The simple formula may not apply in the immediate 
neighbourhood with sufficient accuracy, but equation (13) can be used with success 
if suitable values of 4 can be found. These A values will differ from point to 
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point, being greatest for the two or three squares (or groups of squares) adjacent 
to the singularity, and falling very rapidly as the singularity is left behind. Some- 
times the As can be estimated by comparison with a known solution for boundaries 
which, in the immediate neighbourhood, are similar. Thus for the flow round, 
or in, a corner, in an otherwise infinite field there exists the simple form 


w=kz" 


The value of n is fixed by the angle of the corner. Information regarding a few 
cases is given in Fig. 6 and the values of K given there can be used close to such 
corners, even if the other boundaries are very different. 


The data in Fig. 6 were deduced as follows. Actual values of the function 
considered were calculated at all lattice points near the corner by using the 
appropriate form (Column 4) deduced from (14). It was then found that for a 
point, e.g. B (Case 1, Fig. 6), the “twelve” formula gave a value which had to be 
multiplied by a factor K=1-018 to make it correct. Either this factor can be used 
or from it an appropriate A can be deduced. The value of K is independent of 
the size of the square, but similarity must be maintained in the relation of the size 
of the square to its position. Note that in the case shown the value of ¥ on the 
surface is zero. If, in the actual problem, it is not zero, K is simply applied to the 
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difference between the value of ¥ and its value on the surface. Usually in the 
cases shown the A values are negligible unless the group of squares actually touches 
the corner. Thus only three As are required for the “twelve” or “twenty” 
formulae. If desired, suitable As can easily be deduced for use with single squares 
or diamonds. 


For operating on log q and @ near a stagnation point values of A will be found 
in Ref. 8. 


5.3. Example of the use of As 


It is proposed, as an example, to study the flow of a perfect fluid in a two- 
dimensional parallel-sided channel at a point where there is a 90° bend (see Fig. 7). 
This will be done by three distinct methods. The first method is the most obvious, 
namely operating in the x, y-plane on the stream function y. 


Figure 8 shows the final solution; ¥ is zero on one wall and 1,000 on the other. 
Divide the area into squares as shown. Then the y values to the right are 250, 
500 and 750. Start by writing guessed values at the other points. There is 
symmetry across the diagonal AB, so beyond this line nothing is written. Squares 
and groups of squares have to be straddled over the symmetry line when dealing 
with the points on it and near it, but the values when needed are easily picked out 
from the appropriate image points. 


The fastest method of improving the guessed values is to begin operating by 
using the “hundred” formula. The reason is that while it uses sixteen values at 
each application, five are on the top boundary, fixed at 1,000, and five on the lower, - 
fixed at 0. Thus only six are affected by errors in the guessed values. Having 
worked along the centre line the “twelve” or the “twenty ” formula is used to fill 
in the other two lines. 


Then, starting again on the centre line, the method of differences, already fully 
described, is used to complete the solution. The question will be asked: “ How far 
to the right must the 250, 500 and 750 be placed?” These values can be placed 
initially four or five squares from the corner. In carrying out the operation on the 
differences, these differences are allowed to spread as far as they will. Anyone 
actually doing the process will have no difficulty with this. It is necessary merely 
to keep writing zeros on the lines beyond and the limit establishes itself. 


The values of A appropriate to the “twenty” formula are shown in position 
at the points where they differ significantly from zero. In the later stages of the 
solution the effect of these on the “hundred” formula will need to be taken into 
account. This is done approximately by the reciprocal theorem already enunciated. 
In Fig. 20 the effect has been calculated of A=1 placed at the centre of sixteen 
squares when the “twenty” formula is used. This was done by settling the field 
with zero on all boundaries and A=0 everywhere except at the centre. By the 
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Fig. 8. 
Stream function ¥ at a 90° bend, ¥ on x, y. 


Reciprocal Theorem (Section 5.1) it is known that had A=1 been placed at any 
other point, say A, and zero everywhere else, the effect at the centre would have 
been 0-232. 


This provides a method of summing up the contributions of the As at the 
nine internal points of a square such as that in Fig. 20. The total effect at the 
centre is SFA where the F factors are taken from Fig. 20 and the As are those in 
use for the “twenty” formula—in fact, the values already written in Fig. 8 (—3, 
—13 and —2). Thus suppose the “hundred” formula is applied to the point D. 
There are two A values inside the large square centred on D (Fig. 8), namely — 13 
and —2. Hence the value given by the “hundred” formula (446) is increased by 
which is 


{2x 0-408 + 13 x 0-232} = 


and so 442 is obtained. 


There are only five points where the correction as calculated by this method 
has a value; on writing these in the work can proceed. 


Suppose it is desired to test if remaining errors due to neglected terms have been 
serious. One method is to use a finer grid in the vicinity of the corner. The details 
are shown in Fig. 9, in which the lower line of figures is taken from Fig. 8. Squares 
and diamonds are used alternately for the next line and the “twenty” formula for 
the rest. On working over the field no change has taken place in the three pivotal 
values 460, 613 and 698. Thus there is nothing seriously wrong. Probably if 
another figure were retained some movement would result, but it seems that Fig. 8 
contains values of ¥ which are reasonably accurate. 
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Fig. 9. 
Detail of splitting the grid at A. 


6. The Conjugate Function 


When one of a pair of conjugate functions is known mathematically the other 
can be found. The corresponding statement for the present purpose is that, when 
values at all lattice points in a field of one of a pair of conjugate functions are 
known, the values of the other can also be found. 


Thus given the values of ¥ on an x, y-grid, ¢ can theoretically be determined | 
by using 


B B 


ov ov 
A A 


The differentiations and integrations are to be done numerically. Difficulties how- 
ever very often occur. The most frequent arises from the well-known fact that it is 
impossible to obtain the derivative accurately at the end of a column of figures, 
except in the special case where anti-symmetry can be used and, strictly speaking, 
this is not then the end of the column. For this reason the lines along which the 
integrations are made should avoid, as far as possible, the edges of the field. 


The ordinary formulae of numerical differentiation take no account of the 
fact that the field is orthogonal and so are not as accurate for Laplacean fields as 
they might be. The following method takes full account of the field’s properties 
and so would seem to be the most accurate available. 


Suppose that at a point O in a field in the x, y-plane it is assumed that ¥ in the 
neighbourhood can be expressed by Taylor’s expansion. Differentiate this with 
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respect to x and integrate with respect to y. Making use meanwhile of Laplace and 
its derivatives (Section 4.2) the following simple expression is obtained: — 


+1 
00 


This can be taken as the fundamental expression by means of which the integrals 
can be built up step by step. 


One method of using (16) is as follows (Fig. 10). Consider the values of ¥ at 
A,B,C, spaced at intervals of 2 along the line y=1 and the values of ¥ at A,B, and 
so on along y= —1. 


Express these in Taylor’s expansions from the origin O and solve for the x 
derivatives required in (16). Substituting these in (16) gives the integral along 
A, A, as 

48 /=14(B, + B,— B, ~ 


If the origin is on the field edge a suitable form is 


or more accurately 
241=39(B,+ B,)—21(C,+C.)+5(D,+ D,)—23(A,+A,) (19) 


The first and last of these take full account of third order derivatives. 


x,y plane %.v plane 


Fig. 11(a). Fig. 11(b): 
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6.1. Applications 

The formulae just given could be used to deduce values of the potential 
function from those of the stream function in Fig. 8, but if the ¢ and vy lines are 
required it would still be necessary to cross plot along x and y lines as a preliminary. 
An entirely different method of doing the problem is as follows: — 


If w=f, (z), where w=9+i) and z=x+iy, 
then z=f,(w) and V7x=V?y=0 . . 


where we are now operating in the 9, y-field. 


Suppose y is to be found at every intersection point in Fig. 7. The ¢,¥-field 
is used as a grid on which to operate on y. Fig. 12 shows the final stage of the 
solution. To the right of the symmetry line AB (see Fig. 11(a) and 11(6)) y=0 
is a lower boundary and y=800 is an upper (the channel is assumed arbitrarily to 
be 800 units wide). To the left of AB the boundary values were first obtained from 
a freehand sketch. These were gradually improved from the following considera- 
tions. In Fig. 11 two points p and / are taken, symmetrically placed with regard 
to AB. Then if the origin is at B the x co-ordinate of / is the y co-ordinate of p. 
So if, from the data to the right in Fig. 12(a), the conjugate field is constructed 
(e.g. by using (17), (18) and (19) ) the x co-ordinates to the right are obtained. That 
is, the y co-ordinates to the left are known and so new boundary values can be 
written in to the left. The process can then be repeated, settling the field each time. 


On the right side of Fig. 12(b) are the x increments deduced from the right 
hand side of Fig. 12(a@) by equations (17) and (19). These increments are now 
treated as y increments and used to build up the left side. Thus starting with 626 
at G, adding 161 gives 787 and so on. It will be seen that there is very little change 
from 12(a). The next step would be to resettle the field to the new boundary values, 
using the “hundred” and the “twenty” formulae. As appears from Fig. 6 it is 
necessary in this solution to use As in both the inner and outer corners. 


The whole process is convergent because the boundary values to the left have 
little effect on the field to the right compared with the fixed boundaries on that side. 


Accepting the values shown as sufficiently accurate, the field (Fig. 7) can be 
plotted straight off, because figures to the right are y values and figures to the left 
are the corresponding x values. 


6.2. A third solution 


Since log1/q and @ are conjugate functions, either or both of these can be 
operated on. It so happens that a very rapid solution to the problem just considered 
can be obtained by using @ and working in the ¢, Y-plane. 


The boundary values are definite, as shown along the top and bottom of Fig. 13. 
But only one side, the right, need be considered. Along AB, 6==/4, but when 
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Fig. 13. 


operating squares to the right it makes matters much easier to write =/8 at A and 
B. The reason will be fairly evident if the conditions in the corner of the square 
adjacent to A are examined. 7/8 is in fact the average at the corner of the gradients 
on the two sides of the square. 


On trial it is found that using this artifice the values of A are negligibly small. 
The 6-field is shown settled in Fig. 14. If values of the velocity g are required 
they are found by obtaining the function conjugate to @, that is log1/qg, by the 
formulae discussed in the previous section, care being taken to avoid the corners 
where log 1/q is infinite’. 


If the co-ordinates are required these can then be found by using 


_ [cosé = 
x= | and y= q dt (21) 


It is evident from the symmetry of Fig. 13 or 14 that along the line ~=2, 00/dy 
is zero. Hence 0 (log 1/q)/0¢ is zero, and so q is constant along y=2. It is thus 
particularly easy to integrate along this line and to obtain x and y. The results 
are shown at the bottom of Fig. 14. The integration starts at the right with y= 400. 
At the line of symmetry it is 625 (cf. 626 from Fig. 12) which, from the symmetry 
of the problem, forms the starting point for the x integral shown on the next line. 
The values can be compared with those obtained in Fig. 12(a) by a totally different 
method. The agreement is as good as can be expected, bearing in mind that the 
solution given in Fig. 12 was not complete. 


7. Viscous Flow Equation 
When attempting to deal with equations of the type 


it is advisable to operate simultaneously on two fields, e.g. stream-function and 
vorticity" Write the equation as 


V% =f (x, y) 
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Values of ¢ and ¥ are guessed throughout the region and operated on separately. 
The boundary values of ¥ are presumably known, but the boundary values of ¢ 
must be found as the solution proceeds. 


At a solid boundary the value of ¢, can be found from'2? 


where Y;=the value of ¥ on the boundary and y,,¢, are values of ¥,$ at a point A 
distant n, from the boundary. 


Details for the case of a finite Reynolds number will be found in Refs. 1 and 2, 


but the method of obtaining (, given there is probably inferior to Wood’s method 
just given. 


Here it is only proposed to deal with the case when Reynolds number tends 
to zero. 
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8. VY¥=0 


To keep the formulae identical with those in Refs. 1, 2 and 4, write this 


The working formulae then become? 


“ Twelve”: 


“Hundred ”* : 


In problems of this kind it is nearly always desirable to obtain first a solution of 
V*n= V?£=0 for the given boundaries where €, are the potential and stream 
functions for the irrotational non-viscous flow between the boundaries. The 
£,7-grid then forms a lattice on which the viscous solution can proceed (see Ref. 1 
for proof) without the necessity for troublesome interpolations at the boundaries. 
When this is done the value of n in equations (27) and (28) must be divided by g 
(the velocity in the Laplacean solution), this being a measure of the scale of the 
grid at the point. In the same way the value of n, in (24) must also be divided by 
the local value of gq. 


Problems of this kind usually take a great deal of work. Very often certain 
parts of the field demand a very fine grid, so fine that tens of thousands of squares 
would be needed to cover the whole field. Accordingly the greater part of the 
field is worked with a wide network which is subdivided as the difficult points are 
approached. Often this subdivision is made three, four or even five times, so that 
the smallest square used may be less than one thousandth of the area of the largest. 


It was found difficult to get a problem which would allow of a solution short 
enough for inclusion here, but the following somewhat artificial example brings out 
the essentials of the method without necessitating subdividing. 


8.1. Example 

Consider the slow two-dimensional flow of a viscous fluid between parallel 
planes four units apart. The mean velocity at infinity is taken as 64 units per 
second. The disturbance at a constriction produced by a small symmetrical swelling 
on each wall, as shown in Fig. 15, is studied. 


The first step will be to solve the perfect fluid case to obtain the ¢,-grid on 
which the viscous solution will be made. Here a velocity of unity is assumed in 
the undisturbed channel so that the boundaries are y= +4 and »=0. The length 
of the swelling is about 4 units. The height h is left general but is assumed to be 


*The assumption is made here that the last term has the same value as that deduced in Ref. 4 
for use with (8), and this cannot be seriously wrong. 
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small. The shape of the swelling is defined by specifying the slope of the wall at 
specific values of €. It is true that if A be finite the € lines will be pulled together 
somewhat, with a consequent shortening and steepening of the swelling, but this 
is a second order effect. 


Take the wall slopes to be 
-2 -1 © 4 2 


where @ is the inclination of the surface and y is a small angle. The height of the 
swelling is found by integration to be h=y. There is symmetry about »=2 and 
about €=0 so only one quarter of the whole field need be used. In Fig. 16 the given 
gradients are written as boundaries in the €,7-field and the field is settled. The 
conjugate field (log 1/q)/y is then built up. Fig. 17 shows values of (log q)/y and 
the paths of integration along which they were formed. 


Put L=(logq)/y. Then aoe or, since y is small, 
=1l+yL ‘ ‘ (29) 


Thus the required qg at every point is asi in terms of en small angle y. At 
this stage a definite small value could be assigned to y and (26), (27) and (28) could 
be used to improve guessed values of ¥ and ¢. It is proposed, however, to retain 
y algebraically and to operate on the perturbations brought about by the distortion 
of the walls. 
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029 O09 O05 O02 ro) + Integration along 
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Fig. 17. 
(log g)/y on &, 1. 


At each point replace y by ¥+p’ and ¢ by (+ where vy and ¢ are the undis- 
turbed values. Noting that (26) and (27) apply both to the disturbed and to the 
undisturbed case, putting n for n/q and writing p’ =p, v’ =v, the working formulae 


12 Sve +2 Sv, 

12 po==pc+2 — 32 + 64 n%XoLo (30) 
are obtained where suffix O refers to the centre of the group and C and a refer to 
the corner and mid-side values respectively. 

The end terms have been obtained by putting 
n? 
Similarly the boundary value of v is given by 


where v,,=boundary value of v and suffix 1 refers to values at a distance n, (in the 
£,-plane) from the wall. 


~ n®(1—2yL). 


The fields, completed to two figures, are given in Fig. 18. The values enclosed 
in rectangles are values of (ZL which do not change as the work proceeds. Values 
of p=p’/y are above and values of v=v’/y below the rectangles. The undisturbed 
values of ¥ and ¢ are given at the right. 


Values of ¥’ and ( have now been obtained at all grid points in the form 


V=v+yp and 
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Fig. 18. 
Viscous flow through a contraction (R->0). 
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the dash indicating the perturbed value. 


The loss in pressure is obtained by evaluating along the centre line or along 
the wall the expression 


[55% 


but since q= 3 = ud this can also be written 


iE dé, 


OC+vy) _. (av 
which becomes | dé (2) 


Since V*(=0, equations (17) and (19) are available for evaluating the integral. 


It appears from the details given in Fig. 18 that, for the second term in (32), 
2 x 25y is obtained along the centre line and 2 x 23y along the wall. Taking the 
mean to be 48y it is found that the equivalent length of passage giving the same 
loss as the constriction is 2y, i.e. (y/2)x width. An identical result is obtained by 
neglecting the change of velocity profile and simply assuming the pressure gradient 
to be proportional to (mean velocity)/(width). This in fact constitutes a partial 
check on the whole calculation. 
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THE ARITHMETIC OF FIELD EQUATIONS 


To find the velocity at any point on the plane of symmetry :— 


u mvelocity= = q 


= (5 +752) 


oP ) 
(1 + +yL 
For example on the centre line at the plane of symmetry 
u’ = 96 (1+ 0-34 y). 


As mentioned already, there is nothing to prevent this type of solution from 
being extended to a finite constriction at a finite Reynolds number. The first order 
solution used as an illustration was adopted to save time and space. 


9. Propagation of Error in Laplacean Fields 

It is interesting to examine, by taking specific examples, how an error at one 
point affects the field. This error may be due to a mistake, to neglected terms, or 
to rounding off the last figure. 


Figure 19 shows (F) the effect on a group of sixteen squares of a A of unity at 
the centre square only. By the reciprocal theorem already enunciated the coefficient 
F at, for example, D shows the effect at the centre of an error of unity at D. 


Suppose that a field of this size is being worked using the unit square. Perhaps 
the guessed values have been everywhere low and the field is gradually rising to 
equilibrium. If it is decided to stop when each square is within exactly one unit 
of satisfying (4) then the total error at the centre is obtained by taking the sum of 
all the coefficients on Fig. 19, i.e. %(F), which is found to be 9-68. In Table I 
values of &(F) are given for a number of square fields. s? is the number of unit 
squares, e.g. s*=16 for Fig. 19. It is seen to be near enough to say that =(F) 
is 0-6 

Thus as the size of the field goes up it becomes increasingly important to 
make sure that the process has been carried far enough. If the field being worked 
is not square a rough guide is to draw on it the largest possible square group of 
squares or diamonds, and to take its area (in terms of the side of the unit square) 
as s? in the rule just given. 


In Fig. 20 the same size of field has been treated in a similar manner, but the 
“twenty ” formula used. It appears that =(F) is only about 2/5 of the value found 
in Fig. 19. This is, as already mentioned, one of the advantages of using group 
formulae. 


The foregoing discussion refers to a systematic type of error always of the 
same sign. Now consider the effect of random errors, e.g. those introduced by 
rounding off the last figure. The error introduced into any value ought never to be 
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Fig. 20. 
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greater than 0-5 in units of the last figure. If it so happened that all were of the 
same sign, the cumulative effect would be 4 (F), but as positive errors are as 
likely as negative it appears on statistical grounds that the standard deviation of the 
centre value may be written 


where k has a value not far different from 1/ / (12). 


TABLE I 
s? = (F) 0-6 s? > (F?) 1+ 3s? 


1 

4 2°67 2-222 

9 4:184 

16 9-68 6°866 

25 14-98 10-276 

36 21-46 14-414 

49 29°12 19-287 

64 37:96 24:90 
By considering the known solution for the deformation of a flat square plate loaded 
uniformly it is found that when sco, =(F)>0°588s2. 


Values of = (F?) have been calculated for a number of cases (Table I). The 
last column is 1+ s? and it will be seen that throughout the range considered and 
probably well beyond it this gives a sufficiently good approximation to = (F’). 


Thus, for a large square or diamond the rough rule for the standard deviation 
at the central point due to rounding off errors is: — 


the unit being the last retained decimal. A method of applying this to fields which 
are not square has already been suggested. 


When using group formulae the effect is less and it is again tentatively suggested 
to replace s? by the ratio of the area of the largest square which can be drawn in 
the field to the area of group embraced by the formula being used. 


Suppose it is proposed to work on a field for which s?=324, and the central 
value is required to have a standard deviation of less than 0-001. By (35), c= +3 
in the last place. Thus it ought to be sufficient to retain only one figure more than 
is required, i.e. the solution is continuously refined until it is settled to +0-0001. 
This, as has been explained, is quite different from refining until each value has 
approached to within 0-0001 of the value from the formula being used. In the 
latter case there is a possibility of an error of = (F) x 0-0001 or 0-6 s? x 0-0001 which 
is 0-019 or 19 times the permissible error. 


This example illustrates how extremely important it is in a large field to ensure 
that the values are really settled. It also shows the advantage of, say, the “ hundred ” 
formula over the unit square formula. The former brings down the effective value 
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of s to one quarter of its value for the latter, and so effects a commensurate increase 
in accuracy. 


This is not a hypothetical case. In one actual problem worked by a colleague, 
eight fields each with s?=256 were settled satisfactorily by the equivalent of (28). 
A subsequent examination showed that the actual errors seldom exceeded 2 in the 
last place. 


10. Conclusion 


In conclusion it should be pointed out that the formulae and methods given in 
the paper are intended for desk work. With electronic calculating machinery 
designed for this type of work it is almost certain that the original Liebmann method 
would be the most suitable. Nevertheless the use of the double field for V4v 
equations would still be an advantage. 


It has been asked why formulae for the centre value making use of points 
internal to the large squares are not advocated. The answer is that it seems that 
the larger the clear area covered by each operation the greater the movement 
obtained and so the faster the solution proceeds. Anything inside the square will 
tie down the movement and so reduce the speed. 
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On the Low Aspect Ratio Oscillating 
Rectangular Wing in Supersonic Flow 


JOHN W. MILES 


(University of California) 


Summary: The Laplace transform of the lift distribution on an oscillating 
rectangular wing in a supersonic flow is obtained by separating the linearised equation 
for the velocity potential in elliptic (cylindrical) co-ordinates. The results for the case 
of no spanwise distortion are expanded in ascending powers of the aspect ratio in order 
to compare with the slender body theory, and the longitudinal stability derivatives are 
calculated. It is found that at either supersonic or transonic speeds single-degree-of- 
freedom instability in pitch is impossible insofar as the fourth power of the aspect 
ratio is neglected. 


1. Introduction 


The question of supersonic flow over oscillating wings of vanishingly small 
aspect ratio has been discussed in a previous paper’, the appropriate treatment 
therein being designed as “slender body theory.” In the case of a rectangular 
wing in steady, supersonic flow, Stewartson’’ has succeeded in bridging the gap 
between the results of slender body theory and the well known results for wings of 
large aspect ratio. In the following, an attempt is made to bridge this same gap 
for the oscillating wing. 


The boundary value problem posed by the rectangular aerofoil in supersonic 
flow is analogous to the problem of acoustical diffraction by an infinite, plane strip, 
with the streamwise co-ordinate in the former case corresponding to time in the 
latter. This diffraction problem was first treated by Schwarzschild’, and his pro- 
cedure has been utilised by Gunn? to treat the steady flow over a rectangular wing 
of effective aspect ratio (GA’) greater than one half. Similarly, Stewartson‘’ has 
handled the problem of the oscillating, rectangular wing of effective aspect ratio 
greater than unity by analogy with Sommerfeld’s classical problem of diffraction by 
a half plane. 


For wings of effective aspect ratio less than one half, as indeed for the original 
diffraction problem, Schwarzschild’s technique becomes unwieldy. To circumvent 
the attendant difficulties, Wien‘ and Sieger‘*) proposed the separation of the wave 
equation in the co-ordinates of the elliptic cylinder; this procedure led to useful 
results only much later (with the availability of tabulated Mathieu functions) in the 
work of Morse and Rubenstein’. The introduction of elliptic co-ordinates in the 
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aerofoil problem leads to a solution that is formally complete for any size, but this 
paper is concerned primarily with the wing of small aspect ratio in order to delineate 
more explicitly the proper regime of the slender body theory. 


Notation 


~ 


M 
U 
b 
f 
h 


~ £2 9 ~ 


is) 
= 


Qe 


effective aspect ratio (8A’); also Laplace transform of « 
aspect ratio (see equation (26) ) 

see equation (23) 

Euler’s constant (0-577 

lift coefficient 

moment coefficient 

see equation (22) 

modified, non-periodic solution to Mathieu’s equation; see Ref. 10, 
p. 155 

Laplace transform of h (x) 

Mach number 

free stream velocity 

semi-span of wing 

sonic velocity 

displacement of point on aerofoil 

influence function for potential; see equation (9) 
influence function for spanwise average of potential; see 
equation (11) 

imaginary unit 

reduced frequency 

wing chord 

2/4 

perturbation velocity 

cylindrical polar radius 

normalised Mathieu function; see Ref. 10, p. 24 
time 

cylindrical elliptic co-ordinates 

Cartesian co-ordinates 

effective incidence; see equation (7b) 

(M? — 1} 

see equation (32) 

kM 

b [8?p? + 2ikM?p — k?M?}} 

mass density of air 

dimensionless velocity potential 

spanwise average of ¢ 

angular frequency 

real part of 


A 
A’ j 
A,” 
CL 
Cy 
Gek,, 
A 
u,v 
a 
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2. Formulation of Problem 


Consider a rectangular wing of infinitesimal thickness and semi-infinite chord, 
bounded by x>0 and | y| <b, that is placed in a supersonic flow of velocity U 
and executes a small, monochromatic motion, specified by 


z(x,y. Ale f (x, y)] (1) 


with respect to its equilibrium position, z=0. The right handed, Cartesian co- 

ordinates (x, y, z) are assumed to be dimensionless with respect to a characteristic 

length /. If, in accordance with the usual idealising assumptions, the resulting 

perturbation velocity q is expressed as the gradient of a dimensionless potential, ¢, 
assumed to be everywhere small. 

B=(M?-1} . ‘ (6) 


where k is the reduced frequency and M the Mach number. 


The boundary conditions to be satisfied by ¢ are 


¢:(x,y,0+)= — a(x, y)=f. (x. y)+ikf (x,y), x>Oand|y|<b . (7b) 
o(x,y,0)=0, x>Oand|y|>b . (Tc) 

| (x,rcos @,rsin@)|<K, r—>oo Ce 


Equation (7a) implies that no disturbance can exist forward of the aerofoil and is 
a direct consequence of the hypothesis of supersonic flow. Equation (7b) expresses 
the requirement that the flow be tangential to the aerofoil, the symmetry of this 
boundary condition implying the anti-symmetrical character of ¢ with respect to 
z=0. The latter consideration, together with the requirement of continuity of 
pressure except across the aerofoil, establishes equation (7c). Equation (7d), where 
K is a constant, guarantees the correct behaviour of the solution at large distances 
from the aerofoil. 


The pressure at the upper surface of the wing is given by the linearised 
Bernoulli equation, the result being 


p(x, y, —p.U? {e [92 (x, . (8) 
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7 
U=0 V=0 
Fig. 1. 


Plan and chordwise views of low aspect ratio aerofoil, showing Mach 
lines from leading edge corners and successive reflections. 


where p, is the mass density of the undisturbed air. For the calculation of this 
pressure it is convenient to seek a solution in the form 
b 


—b 


In the important special case where the displacement of the wing is independent 
of y, and where the spanwise pressure distribution is of no direct interest, it is 
expedient to deal with the spanwise average 

b 


1 
—b 
which is to be expressed in the form 


xz 


va= 


0 
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The problem is now to obtain a solution to equation (3), subject to the 
boundary conditions (7). 


3. General Solution 

The difficulties introduced as a result of the specification of the boundary con- 
ditions (7a, b, c) over three distinct regions of the plane z=0 may be effectively 
dealt with by introducing a Laplace transformation with respect to x and posing the 
co-ordinates of the elliptic cylinder defined by 


y=bcoshucosv ; ‘ . (12a) 
z=bsinh u sin v ‘ . (12d) 


In these co-ordinates, the regions y<—b, | y| <b and y>b correspond to v=7, 
u=0 and v=0, respectively. Furthermore, u and v approach log.(r/b) and 4, 
respectively, for large r. Then, denoting the transformed functions by the corres- 
ponding upper case letters, we write 


(p,u, v)= L {¢ (x, b cosh u cos v, b sinh u sin v} = o(x,y,z)dx (13) 
0 


A(p,v)= {«(x, b cos v)} ? ; . (14 

P,,,+P,,=A? (cosh? u—cos?v)P . : 

A= Bb [(p + + = Bb [p + ix (M + [p+ ix (M — 1} . (16) 

®, (p,0,v)= — bA (p, v) sinv (18a) 

(p,u,0 or =)=0 ‘ . (185) 


The branch cuts for A are to be taken from p= —ix(M+1) to p=O—ix(M+1) 
and the phase chosen to ensure % (A)>0 when & (p)>0. In the end results, 


however, these branch cuts affect only logarithmic terms in A. 


A solution to equation (4) satisfying the boundary conditions (7) is given by 


2b Gek, (u, —q) sen (v, — 9) A (p, bsinv)se, (wv, —g)sinvdv (19) 


n=1 Gek,, (0, — q) 


® (p, v)= 


where the notation for the solutions to Mathieu’s equation is that of McLachlan". 
Taking the Laplace transform of equation (9), utilising the Faltung theorem, 
and comparing to equation (19), gives 


G (p, 2, w)=2 F.Q)se, (v, —q)se,(w, —q) . 2h) 
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_ _ Gek, —q) 


The special case where a is independent of y, and therefore of v, requires only 
the odd Mathieu functions, which have the expansion (see Ref. 10, Section 2.18, 
equation (4)). 


Substituting in equations (19) and (10) and comparing to equation (11), gives 
WV (p)=H (p) A (p) (24) 
which completes the formal solution to the problem. 


4. Approximate Representations 


In principle, the results of Section 3 are applicable to rectangular wings of 
arbitrary aspect ratio, but the exact inversion of H(p) in terms of (presently) 
tabulated functions appears as a formidable, if not impossible, task. However, since 
the objective is the study of the low aspect ratio wing, it will suffice to expand in 
powers of 8b (except near M=1, where the expansion is in powers of kb; see later). 


If the characteristic length / is chosen as the wing chord, the parameter Gb is 
directly related to the “effective” aspect ratio A by 


where A’ is the usual, geometric aspect ratio. 


As a preliminary to the desired expansion, the Bessel function representation of 
the Gek, (Ref. 10, Section 13.30, equation (7)) is introduced, the substitution of 
which in equation (22) yields, after some manipulation, 


Fong (A)= (x) { (ar+ 1)+ (5) (3) K, (5) 


where 7 and K denote modified Bessel functions. 


The principal difficulty in effecting an expansion in powers of A lies in the fact 
that each of the F2,4; has at least one pole. These poles will lead to terms of the 
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form exp(—constant x A-'). Such terms are negligible compared to terms like A" 
for sufficiently small A, but Stewartson’s results make it evident that it is essential 
to consider at least the smallest (in magnitude) of the poles if acceptable results 
are to be obtained. Upon expanding F, (A) and retaining terms up to A‘, they are 
found to be located at A=A, and A,, where 


A, = . (28) 


in agreement with Stewartson’s result. (According to Stewartson, the remaining 
zeros lie outside of | A |=4, and the present results tend to confirm this estimate.) 


To complete the development of H, the singularities at A, and X, are separated 
out, namely a_,(A—A,)~', and a_,(A—A,)~', where a_, is the residue. 


a_,=lim (A-—A,)H (p)=0-769 ‘ (29) 
‘>A, 


The remainder, namely H (p)—a_, (A—A,)~*—a_, A —i,)-}, then may be regarded 
as regular in A and expanded in powers thereof. The end result, arrived at after a 
tedious but straightforward development of equations (25) and (27) with the aid 
of the various formulae in Ref. 10 may be cast in the form 


H (p)=1+ *) (3) + x) (3) +O 


"+3 (fama) "(BY 


v=C+log.(*), C=0577... GD 


where C is Euler’s constant, and the phase of y is chosen such that y=A, at «=0. 
Terms involving A*°, although used in the determination of A, and a_,, have now 
been dropped. The terms of the geometric series in equation (31) serve to cancel 
the lower limits (x=0) arising after integration of the exponential terms obtained 
from the inverse transforms of [p+ikM—2A~-' (y, 

The approximate expansion (30) rests on the assumption that A is small but 
otherwise places no restrictions on k and M beyond those already implicit in the 
linearised theory. However, if kKM«1, as it often will be in dynamic stability 
problems, additional simplification may be possible. Three cases present themselves; 
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as (M—1)-'k is small, approximately equal to, or large, compared with unity. In 
the first of these it is convenient to expand H about k=0 and write 


H (p)=H, (p) + ikH, log. (Ap) G3) 
se 


(34D) 


(35) 


In the intermediate case where (M-—1)-'k is near unity no appreciable 
simplification of H results, since there remain terms containing log. (p + constant), 
the inversion of which lead to sine and cosine (or exponential) integrals. However, 
if (M—1)«k« 1, as will be the case for M sufficiently close to unity, 


putting M=1 in view of the approximation already involved. The resulting expan- 
sion of H (p) then yields 


H (p)=1+ 7 loz. . +it | (74) + 


where the poles have been neglected, since they now are located at a distance of 
order 2\,?/(kKA”) from the origin. Separating this as in (33) and neglecting the 
term (— kA”p/32) in H, as negligible compared to unity, gives 


H,(p)=1 
H, (p)= [c+5 
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5. Lift and Moment Coefficients 


The lift and pitching moment (acting about an axis at x=a; positive when 
tending to raise the leading edge) on the wing are given by 


1 bl 


1 bl 


where the wing chord has been chosen as the characteristic length. As measures 
of these quantities, the complex lift and moment coefficients are introduced as 


Calculating the spanwise integral of the pressure at the upper surface (equal and 
opposite to that at the lower surface) with the aid of equations (8) and (10), gives 


CL =4 (ve + iky) . (43) 


Cu=4 | (a—x) +iky) dx. ‘ (44) 


To find the Laplace transforms of equations (43) and (44) the upper limits of 
the integrals are taken as variable, and substituting Y from equation (24) gives 


6. Longitudinal Stability Derivatives 


The results of the preceding sections are now applied to the calculation of the 
longitudinal stability derivatives associated with the motions 
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The Laplace transforms corresponding to unit values of the respective amplitudes 
w and @ are given by 


=p" . P ‘ ; (49) 
Retaining only terms of O(k), the stability derivatives may be put in the form 


with similar expressions for the moment coefficients, after replacing the suffix “L” 
by “M.”) Cyo+ikCzsa at is the lift coefficient resulting from the downwash distribu- 
tion z= 1, neglecting terms of O(k?), while C,, results from the distribution (x—a), 
neglecting terms of O(k). Substituting (49) and (50) in (45) and (46) and comparing 
with (51), gives 


C. tA -2 -1 
Lan . (54) 


{p-* H, (p)+p~? A, (p)} 


where each of the inverse transforms is to be evaluated at x=1. 


The required inverse transforms are 


L (4) +12 [3 + 


+2R | 


4] 
240 


| 
f 

| 
| 

| 


LOW ASPECT RATIO OSCILLATING SUPERSONIC WING 


L -*{p-*H, - [3 +108. (4 Nay (5) (4) + 


L (4) {2 (4) +163 -108. (4) ] (4) + 


L -{p-*H, +1og.(4)](4) + 
+AK| 


The obtaining of further information demands the substitution of typical 
numerical values, some examples of which are given in the following section. 


7. Lift Curve Slope 


To estimate the accuracy of the foregoing approximations, consider first the 
lift curve slope (cf. Ref. 11) which, as given by equations (53) and (59), may be 
written 


loge (4)]4 L-172e one +0: 275 |} 
(64) 


where C;,, denotes the two-dimensional lift curve slope, namely 


(65) 


Cort 
A’-300 B 


The result given by equation (64) is plotted in Fig. 2. 
The exact (linearised) result for A >4 is given by"? 


Cra 


=1-QA)"?, A>l 
(2A) 


4+ (1-4 
2-A\,_ 1+(1-A?}! 
(668) 


At A=1 equations (64) and (66) both give 0-500 (slide rule accuracy), and at 
A=4 the numerical results are similarly indistinguishable. At A=4, the slender 
body theory alone (C,./Ci.,=7A/8, as first obtained by Jones"*’) is accurate to 
better than 2 per cent., but this agreement is to a slight extent deceptive, since the 
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result given by equation (64) continues to oscillate about the slender body result in 
A<}, as shown in Fig. 2. A similar comparison for Cy. shows the centre of 
pressure calculated from the foregoing to be in error by 3 per cent. at A=1. 

Cy. was selected as a basis of comparison both because it is the simplest of the 
derivatives and because of the availability of the results of Refs. 4 and 12. It is 
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Fig. 2. 
Lift curve slope as given by equation (64) compared with the result of R. T. Jones. 


also the least critical, but it is reasonable to expect.that the remaining derivatives 
would show errors of the same order of magnitude as Cy.. Accordingly, the results 
of Section 6 should suffice for the calculation of the dynamic stability derivatives 
of a rectangular wing subject to the limitations 


A<l, k(M-1)"'«1 . . (67) 


Moreover, the latter restriction may be removed by introducing in place of the 
approximate expansion of (33) and (34) that of (30) or (37). The latter case is 
now considered. 


8. Sonic Stability Derivatives 
Substituting H, and H, from equation (38) in equations (53) to (58) yields the 
sonic derivatives 
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Cma=[1- (4) Jen . 


The approximation involved in these results is indicated by the error term in 
equation (37). In this connection, the linearised theory is valid at M=1 if either 
k-? 8% or A’d} is sufficiently small compared to unity, as may be shown by an 
extension of the two-dimensional analysis of Lin, Reissner and Tsien"*?. 


9. Transonic Pitch Damping 


An important question in either dynamic stability or flutter analysis in the 
transonic regions is the sign of the total pitch damping derivative given by 


At M=1, the substitution from equations (72) and (73) yields 


It is evident that the sign of Cy «,a: is negative definite for values of A’ within 
the limits of approximation, although the damping may be quite small for 2 axes near 
the trailing edge (a= 1). 


A similar result is obtained using the results of Section 6 for k(M—1)-' «1. 
This agrees with an analysis"* valid only for A’ > 1, in which it was concluded that 
negative damping (positive Cy a,a:) could occur only for A’ > 1-93. 


It is of interest to contrast the foregoing results with their two-dimensional 
counterparts. The boundary value for the Laplace transforms (see equations (13) 
to (18)) reduces to 

®,,=(2ikp —k?) ® ‘ . (76) 


The solution to (76) and (77) is given by 
z>0 . . (78) 
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where the branch cut for » passes from —4ik to —4ik— 00, and the phase is chosen 
in such a way as to ensure a positive real part. 


Putting z=0 in (78), the required inversion is given by 


The leading term (for small k) in the pitch damping coefficient (based on a unit span) 
is then found to be 


whence negative (unstable) damping in pitch is indicated for all axes forward of 
the one-third chord point. In connection with this last result, it is not legitimate to 
take the limit k=0, since the linearised theory is not valid there. 
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A Generalised Approach to the Local 
Instability of Certain Thin-Walled Struts 


A. H. CHILVER, Ph.D. 


(Department of Civil Engineering, University of Bristol) 


Summary: The problem discussed is that of the elastic local instability of a uniformly 
compressed thin-walled strut composed of a number of flat component plates. The 
strut is essentially “ open ” in cross-sectional form, while reinforcing flanges are attached 
to the extreme edges. An overall stability equation, based on the small deflection theory 
of plate bending, is derived from the conditions which must hold at the common and 
extreme longitudinal edges of the strut. 


It is shown that the critical compression stress induces a mode of buckling which 
involves a whole number of sinusoidal half-waves in the longitudinal direction. Further- 
more, the number of half-waves is common to all component plates of the strut. The 
general stability equation—a zero determinant of high order—lends itself to expansion 
in terms of 4th order minors, which may again be expressed in terms of seven basic 
functions. A knowledge of these functions is sufficient for the solution of any problem, 
however complex, within the scope of the general analysis. 


Application of the basic functions to the solution of three different problems yields 
results which indicate the need for considerable care in the design of thin-walled struts 
with reinforcing flanges. In struts of this type a longer wavelength of buckling is possible 
than is commonly associated with local instability problems. 


1. Introduction 


Local instability is a form of buckling which may occur in compressed thin 
plates and compressed thin-walled struts. During the past twenty years or so both 
aeronautical and civil engineers have realised that the careful use of thin-walled 
members leads to structures which are not only economical in weight of material, 
but also possess degrees of strength and stiffness comparable to those of the heavier 
types of construction. For practical reasons many thin-walled struts consist of an 
assembly of thin flat plates, connected to form a suitable cross-sectional shape. 
In the early days the design of such sections was based largely on the method given 
by Bryan"? for single flat plates in compression. It was soon realised, however, 
that a treatment of the problem in this way is inadequate, for no account is taken 
of the important interaction effects between one component plate and another. A 
comprehensive approach to some simple problems was made by Lundquist®’, and 
Stowell and Lundquist"); the more recent work of these writers has been based on 
the use of the “moment distribution” method developed by Lundquist”. As a 
result of these researches considerable information is already available on the local 
instability strength of the simpler types of struts. Experimental work by Heimer]? 
on light alloy struts and by the author” on mild steel struts shows that the 
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theoretical critical stresses are in good agreement with those found in tests. The 
maximum strength of a thin-walled strut is reached at a stage subsequent to the 
initial local instability of the section; it has been demonstrated by tests) that the 
maximum attainable stress is dependent upon 


(i) the elastic instability stress, and 
(ii) the “yield” or “proof” stress of the material. 


These considerations have led the author to think in terms of a generalised 
method of deriving the theoretical local instability stresses of thin-walled struts. 
The modern trend in thin-walled construction is to use struts in which the 
unsupported or free edges of component plates are reinforced with “bulbs” or 
“lips.” These reinforcing flanges are added in an attempt to increase local 
instability strength. An essential feature then of a generalised approach to the 
problem must be that of dealing adequately with the flexural and torsional effects 
of reinforcing flanges. 


The paper is devoted to the development of such a generalised method of 
analysis. Application of this method to any particular problem results in the 
derivation of elastic instability stresses which are exact, since no recourse is made 
to the approximations common in the “energy method.” The generalised method 
may be applied to struts having unsymmetrically shaped cross sections, since no 


preference is given in the analysis to particular modes of instability. 


Notation 


number of component flat plates 

length of the strut 

applied and critical compression stresses, respectively 
rectangular co-ordinate axes 

breadth and thickness, respectively, of a component plate 
flexural stiffness of a component plate 

modulus of elasticity of the strut material 

Poisson’s ratio of the strut material 

deflection of a component plate 

Fourier coefficient—a function of y only 

an arbitrary constant 

abbreviations for cosh and sinh, respectively 
abbreviations for cos and sin, respectively 

instability functions 

flexural and torsional stiffness, respectively, of a reinforcing flange 
cross-sectional area of a reinforcing flange 

parameters, defined by equations (7) 

parameters, defined by equations (6) and (8) 

parameters, defined by equations (9) 
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TK 


Fig. 1. Fig. 2. 
The general structural member. Unstable form of the strut. 


a,b,c,d,e, | abbreviated constants appearing in the stability equation (28), and 
f.2.h,i,j J defined by equations (19)-(26) 
K determinantal coefficient 
m_ number of half-waves in the m™ mode of buckling 
R,F,B,B’,S modes of edge support, defined in Section 6 
fb, fs, b’s, b’b, em functions, defined by equations (30)-(36) 
ss, sb, bb 
The suffix m refers to the m™ mode of buckling, while the suffix n refers to the 
n™ component plate. 


2. Statement of the Problem 


The uniform compression of a thin-walled strut presents a typical problem in 
elastic instability. The general structural member, Fig. 1, consists of a number of 
component flat plates, the longitudinal edges of which are either common to each 
other or supported by reinforcing flanges. The strut is essentially short, so that 
overall column instability, in one form or another, is not a primary mode of 
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buckling. The unstable form of the strut under end compression involves 
distortions of the type shown in Fig. 2. Whereas the common edges of the com- 
ponent plates remain straight, the reinforcing flanges may deflect with the component 
plates to which they are attached. 


The discussion is restricted to struts of “open” cross section in order to 
simplify the forms of the stability equations. The term “open” implies that the 
wall of the strut is “continuous” or “unbranched,” and that at a common edge 
there are no more than two component flat plates. Other cross-sectional forms 
may be studied by a modified generalised method, but these cases are left for a 
later discussion. 


The general structural member, Fig. 1, consists of N flat rectangular plates. 
The extreme component plates, the Ist and N“, are restrained both flexurally and 
torsionally at the remote edges by reinforcing flanges. The length L of the strut 
is subjected to a uniform end compression stress o, acting both on the component 
plates and on the reinforcing flanges. Co-ordinate axes, Ox,, Oy,, are taken in the 
plane of the n™ component plate; in the buckled mode of the strut there is no out- 
ward bending of the component plates along the common edges Ox,, Ox,, ... Oxw 
(Fig. 2); b, is the breadth and ¢, is the thickness of the n component plate (Fig. 3). 
All component plates are simply-supported along the edges x=0 and x=L. 


Zn 


tn 


Fig. 3. 
Dimensions of the n*® component plate. 


One 


3. The Unstable Form of an Isolated Component Plate 


Along the longitudinal edges of a component plate the boundary conditions 
depend entirely upon the nature of the surrounding component plates. It is essential 
to discuss the stability of an isolated component plate in terms of arbitrary 
boundary conditions along the longitudinal edges. The boundary conditions are 
defined more carefully at a later stage when relating the equilibrium of neighbouring 
component plates and when specifying exactly the conditions at the reinforced 
longitudinal edges. 


In the buckled condition the deflection at any point of the n™ plate is w,, 
perpendicular to the xy,-plane (Fig. 1). For the buckled mode to be a condition 
of equilibrium, 
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in which D,, the flexural stiffness of the n™ plate, is defined as 


where E is the modulus of elasticity and v is Poisson’s ratio for the material of the 
strut. The equilibrium condition holds only if the deflection w, is small in relation 
to the thickness ¢,. The general solution of equation (1) is simplified appreciably 
by taking w, in the form 


[w= Fasin™™ 


in which F,,,, is a function of y, only. This equation satisfies the boundary 
conditions along the simply-supported edges x=0 and x=L. Substitution of 
equation (3) into equation (1) gives 


[F=A,chay+A,shay+A, cn By+A, sn BY]m, ns 


The expression for Fn, is true only if 8n,,>0. The case for which Bn,»<O0 is 
rather particular and requires a special treatment; however, this case is unimportant 
and is not discussed here. The deflected form of the n™ component plate may then 
be written 


w= (A, ch ay+ A, sh ay +A, cn By+A, sn By) sin 6 


4. The Boundary Conditions at the Longitudinal Edges 


The 1st component plate is reinforced at the edge y, =0 with a fx. _e for which 
B, is the relevant flexural stiffness*, C, is the torsional stiffness, and a, is the cross- 
sectional area. Then a consideration of the bending and torsional actions at the 
junction of the lst component i and its reinforcing flange yields the equations 


[ 
Boat | 


for y,=0. At the junction of the N™ eee plate with its reinforcing flange, 
the corresponding viseiggge’ conditions are 


B 
+a 
*This corresponds to the bending stiffness of the reinforcing flange about a centroidal axis 
parallel to the axis Oy,. 
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for yy=by, in which By, Cy and ay are the relevant flexural stiffness, torsional 
stiffness and cross-sectional area, respectively, of the reinforcing flange attached to 
the component plate. 


At a typical common edge, such as O,x, (n~ 1), the conditions which must be 
satisfied are 


(i) there is no deflection in the (n—1)" component plate, giving w,-,=0 
for 


(ii) the angle between the (n— 1)" and n™ plates in the wy,_,- and wy,-planes 
for the unbuckled configuration is maintained during buckling, giving 


ow] _ [ow] 
40 dy J, =0, 


for and y,=0; 


in the buckled condition there is equilibrium of bending actions in the 
Wyn-,- and wy,-planes, giving 


for y,-,=b,-, and y,=0; and 


(iv) there is no outward deflection in the n component plate, giving w,=0 
for y,=0. 


On substituting for w in these eight boundary conditions, the resulting equations 
are most conveniently expressed in terms of non-dimensional parameters, which 


are defined as 


[{2 ] [{% = 


in which n=0 is some convenient reference datum. Evidently 


[m= 
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Then at the extreme edge of the Ist component plate, defined by y, =0, the boundary 
conditions yield the equations 

F . il) 

[A,r? — A, pe — A, — . 

for m=1 to 00. At the junction of the (n—1)" and n™ component plates, the 
boundary conditions yield the equations 

[A (A,p sh p+ A,pch p— sn CN [A (Ap + m,n =, (14) 

(A,r’ ch p+ A,r’ sh p— A,s° cn g— A,S° g))m,n—1 — (Ayr? — (15) 

[A,+Ag]m,n=0, ; . (16) 


for m=1 to 0. At the extreme edge of the N“ component plate, defined by 
yv= by, the boundary conditions yield finally the equations 


[A, (Ech p— ps? sh p)+ A, (€ sh p— ps ch p)+ 
+A, qr cng)m,~=9, . (17) 
[A, (pe sh p+r’ ch p)+ A, (pe ch p+r’ sh p)— 


A, (qesn (18) 


for m=1 to co. To make the handling of equations (11)-(18) less cumbersome they 
are rewritten in the abbreviated forms 


(Ac)m,n-1=0, 


in which the summations are made from 1 to 4 in each case. The equalities between 
the two sets of eight equations are immediately derived by inspection. 


5. The Stability Equation for a Strut of N Component Plates 
For the whole strut 4N equations are derived from the boundary conditions 
discussed in the preceding section; suppose these equations are assembled and 
written in the form | 
[K]m {A}m=0, (for m=1 to 0), . 
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in which [K], is a square matrix of order 4N x 4N, formed from the coefficients of 
(A, A, Az; A4)m,n etc., While {A},, is a column matrix of order 4N and of the type 
{A,,,..-,Aj,2.-.,€tC.}m. The values m=1, m=2,... etc., in equation (27) corres- 
pond to the Ist, 2nd, ... etc., modes of buckling of the strut. As the end compression 
stress o is increased from zero there is a range of stresses over which equation (27) 
is satisfied only because {A},={A},=...=0. This implies that there is no 
instability of the strut in the forms envisaged. But as the end compression stress 
is further increased, o will ultimately attain some value, o,, say, for which all the 
equations (27) are again satisfied; but whereas in some of these equations the 
matrices {A} are zero, in others the determinants | K| are zero. Generally the 
determinants |K| are not simultaneously zero; if, as the compression stress is 
increased from zero, | K |,, is the first determinant to take a zero value, then the 
m™ mode of buckling defines completely the unstable form of the strut. This follows 
from the condition that if | K |,=0, then {A}, may assume a non-zero value. 


It is interesting to note that there is a mode of buckling common to all the 
component plates of the strut, all plates being deflected into the same number of 
longitudinal half-waves. This conclusion has been widely confirmed by 
experimental work. 


For the m™ mode of buckling the stability equation, expressed in terms of the 


abbreviated notation, becomes 


Az) Az, 
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dy, dz, dy, dy, C22 
fii for for far S12 822 
Ago 
Cio Coo Cao 
dio doz Aye C13 Cos Cas 
fie foo faz 813 82s Bas 
Ciz Cog Cas 


Cin Con 


83n 84n 
hewn hyn 


iw isn ign 


jin jon isw jan 


| 
4 
| 
| 

— 

252 


LOCAL INSTABILITY OF THIN-WALLED STRUTS 


6. Expansion of the Stability Equation in Terms of 
4th Order Minors 


The solution of the general stability equation (28) is greatly simplified after the 
determinant is expressed in terms of 4th order minors. Consider the 4th order 
minors defined by 


RB=|a be ad|, RS=|a bc 
SR=\g jl, BR=|e h i jj, 
SB=|g hc dl, SS=|g hc fl, 
BB=|e hc dl, BS=|e he fi, 


In defining these minors the more important rows of each determinant are 
represented symbolically by the single letters R, B and S. Actually R represents 
a restrained edge, B a built-in edge, and S a simply-supported edge. If a minor is 
equated to zero the resulting equation defines the stability of a single plate with the 
appropriate boundary conditions. For example, RB=0 is the stability equation for 
a single plate restrained at the near edge and built-in at the remote edge. On this 
basis it is reasonable to expect that a minor such as RB is equal to its conjugate 
BR in all respects except perhaps in sign. This is found to be the case, and in fact 


RB=-BR, SR=+RS, BS=-—SB. 


On taking into account these equalities the distinct minors reduce to RB (or— BR), 
RS (or+SR), SB (or — BS), SS and BB. 


The general stability equation (28) may now be written in terms of the symbolic 
notations; the expansion of the equation is effected by using Laplace’s method. 
If N is even, the stability equation becomes 


RB, BB, BS, 0 0 

0 BB, SB, BB, BS, - 

0 BS, SS, SB, SS, - 

0 O BB, SB, 
0 
BB... 
BSy-, SRy 


If N is odd the “ upper” terms in the determinant are unchanged, but the “lower” 
terms take the form 


SBy_, BBy-, BSy-, 
SSy-2 SBy-, 
0 SRy 


=0. 
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In every case the stability equation expressed in terms of minors is an N“™ order 
zero determinant. 


7. The Nature of the Extreme Edge Restraints 

The generalised treatment of the problem has resulted in the introduction of 
minors of the type RB and RS. These determinants may be expressed more simply 
in the following way: suppose that as £—>0 and «—> 0 the limiting value of RB 
becomes FB. In the same way RS, SR and BR take the limiting values FS, SF and 
BF, respectively. Suppose, furthermore, that as £—>0O and «—> 00 the limiting 
value of RB/< becomes B’B. In the same way B’S, SB’ and BB’ are the limiting 
values of RS/e,SR/« and BR/<, respectively. In these limiting values of the minors 
F denotes a “free” edge and B’ a “ partially” built-in edge; of these eight new 
determinants only four are distinct, for 


SF=+FS, BF=-FB, BB =-—B'B. 


It may be shown that the minors and their limiting values are related by the 
equations 
SB BB —BB 
SS BS — BS (29) 
SS SB SB 
BS BB BB 


All these equations show consistency in symbolic representation, since the relations 
SR=+RS and BR=—RB may be deduced easily. 


8. The Basic Functions 
The general stability equation has been expressed firstly in terms of the 4th 
order minors of the type RB, RS, SB, SS and BB. It has been shown subsequently 
that RB and RS may themselves be expanded in terms of FB, FS, SB, SS, BB, B’S 
and B’B. If these latter determinants are written independently of A and yp, then 
seven basic functions are formed. A complete knowledge of these functions makes 
possible the solution of any problem. The basic functions are defined as 
fb= FB/A=(p’s* — q’r*) sh p sn gq — pq [2r’s? + (r*+s*)chpeng], . (30) 
bb= BB/\*? =2pq (1—ch pcn q)+(p? — q”)sh psnq, G4) 
b’b= — pq(r? +s*)(psh peng+qchpsnq). . (6) 
A few numerical values of the basic functions are given in Table I; the values have 
been selected to indicate the general forms of the basic functions. 
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Fig. 4. Fig. 5. 
Angle sections (N=2). Z- and channel-sections (N =3). 


9. Some Particular Cases of the Stability Equation 

The particular cases discussed are those corresponding to N=2, 3, 4 and 5. 
When N =2, the type of cross section implied is commonly called an angle section 
(Fig. 4). The stability equation becomes 


BR, 
| RS, SR, 
If both component plates are geometrically equal, the relevant solution of the 


equation is RS=0; in this case the component plates are in effect simply-supported 
at the common edge. 


=0. 


For the case in which N=3 there is a wide range of cross-sectional forms. The 


more common shapes are known as Z- and channel-sections (Fig. 5). The general 
stability equation reduces to 


RB, BB, BS, 
x 
| 0 BR, 


When N =4 the cross-sectional form is that of a reinforced angle strut (Fig. 6) 
in which the reinforcing flanges are sufficiently thin to be considered component 
plates. The stability equation is 


RB, BB, BS, 
RS, SB, 
0 BB, SB, 
-0 BS, SS, SR, | 
For the case N=5 the cross-sectional form is that of a channel- or a Z-strut 
with reinforcing flanges which again form component plates in themselves (Fig. 7). 
These members are sometimes referred to as “top-hat” sections. The stability 
equation is 


BB, 


| 
0 0 | 
0 0 | 
BB, SB, BB, BS, | =0. 
0 BS, SS, SB, SS,| 
0 0 0 BR, SR,| 
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Fig. 6. Fig. 7. 
Reinforced angle section (N =4). Top-hat section (N=5). 


10. The Complete Solutions of Three Particular Problems 


To demonstrate the use of the generalised method, three typical structural 
problems are solved with the help of the basic functions. In the first example the 
cross-sectional form of the strut corresponds to a “square” symmetrical channel 
of uniform wall thickness (Fig. 8). The Ist and 3rd component plates are quite 
free at the remote edges. If the component plates 1, 2 and 3 are equal in all respects 
then, in terms of the basic functions, the stability equation becomes 


in which the relevant solution is the second root. Values of (c.,b?t/D)! and L/b 
satisfying equation (37) are given in Fig. 8. 


The second example is taken from the first by adding reinforcing flanges to the 
extreme edges of the strut. The reinforcing flanges are narrow plates of the same 
thickness as the walls of the strut (Fig. 9). If complete symmetry of the cross section 
is preserved, then the stability equation is 


RB’. SS+2RB. RS. SB— BB. RS?=0. 


The torsional stiffnesses of the reinforcing flanges are negligible, so that RB and 
RS may be taken with sufficient accuracy in the forms 


RB=FB+(é/)SB, RS=FS+(é/,)SS. 


In terms of the basic functions, the stability equation then becomes 


(38) 


in which the relevant solution corresponds to the second root of the equation. The 
critical stresses for a range of different reinforcing flanges are given in Fig. 9 for 
b/t=100. 
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4 


Fig. 8. 
The buckling of a “ square ” channel strut. 


b, />> dour 0-20 


6 L hb 10 
Fig. 9. 
Buckling of a square reinforced channel strut for which b/t= 100. 
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2 4 5 
L/b, 
Fig. 10. 
Buckling of a square channel strut with wide reinforcing flanges. 


In the third example the reinforcing flanges are sufficiently wide to form 
component plates within themselves. If, in the general strut, plates 1 and 5 are 
equal and plates 2, 3 and 4 are equal, then the stability equation becomes 


+ 1-355 BB/SB (39) 


SS \3—SS BB/SB? 


This equation is true for a strut of uniform wall thickness; the second root is the 
relevant solution. Critical stresses for a range of values of b,/b, are given in Fig. 10. 


11. Conclusions 


It is evident that all problems of the general type are soluble in an exact and 
elegant manner. The complexity of any particular solution depends largely on the 
nature of the cross-sectional form of the strut. 


In general there is found to be a common wavelength of buckling for the com- 
ponent plates in the longitudinal direction of the strut. The number of longitudinal 
sinusoidal half-waves is the same for. each component plate. 


For struts without reinforcing flanges a minimum elastic instability stress is 
reached after a certain length of section is exceeded. This length is usually so short 
for practical struts that the minimum value of the elastic instability stress may be 
considered as the characteristic buckling stress. For struts with reinforcing flanges 
care must be taken in the selection of flange dimensions to ensure that the reinforcing 
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flange is not a major source of weakness. Since struts of this latter type are 
becoming widely used in structural engineering, it is evident that more information 
is required on the behaviour of reinforcing flanges. 
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The Distribution of Lift over the Surface of 
Swept Wings 


D. KUCHEMANN, A.F.R.Ae.S. 
(Royal Aircraft Establishment) 


Summary: A simple method is given for calculating, in incompressible flow, the 
spanwise and chordwise loadings on thin swept wings with symmetrical aerofoil sections 
and of moderate and large aspect ratios. The present note is based on, and continues, 
the investigation described in an earlier paper. The design of swept wings with given 
spanwise and chordwise load distributions, obtained by modifying the plan form of the 
wing or by adding camber and twist, is briefly considered, as well as the design of swept 
wings which give minimum induced drag. 


1. Introduction 


The distribution of lift over the surface of a given wing may be determined 
theoretically by replacing the wing and its wake by a system of vortices. The lift 
force at a point on the wing surface is directly related to the direction and strength 
of the vorticity vector there, and this vector can be found from the boundary 
condition that the velocity component normal to the wing surface at any point, 
induced by all the vortices is equal and opposite to the normal component of the 
main stream. 


On a straight or sheared wing of infinite span, the vortex lines are straight 
and parallel to the leading edge and, since there are no trailing vortices, they alone 
must produce this downwash. The lift is uniform along the span, i.e. the same lift - 
is produced at a given angle of incidence at each spanwise station (corresponding 
to a certain value of the “lift slope,” C,,/a), and the chordwise distribution of the 
lift is the well-known “ flat-plate” distribution for a thin symmetrical aerofoil. 


The position is changed radically if the wing is swept (with a kink at the line 
of symmetry). Experimental evidence, such as that shown in Fig. 1 for a 45° 
swept-back wing of infinite aspect ratio, indicates a considerable reduction of the 
lift in the central region. There are now trailing vortices in the wake (Fig. 2) but 
these produce an upwash and are, therefore, not responsible for the reduction in 
lift. This can only be explained by assuming that the spanwise vortices (i.e. the 
component of vorticity parallel to the leading edge) produce a larger downwash at 
the centre than farther out on the “sheared part” of the wing; and this in turn 
implies that the sectional lift slope varies along the span of the wing and is smallest 
at the centre section. Further, there is a change in the variation of the downwash 
along the chord, which in turn leads to a re-distribution along the chord of the 
spanwise vortices and hence of the load, and this is, in fact, borne out by experiment 
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(Fig. 3). The chordwise loadings at various stations differ in a characteristic manner, 
with a tendency for the aerodynamic centre to move backwards as the centre line 
is approached (Fig. 1). This implies that the vorticity vector curves round in the 
central region, as indicated in Fig. 2. 


On a wing of finite aspect ratio, there are streamwise or trailing vortices 
spreading inwards from the tips. They produce a downwash which reduces the 
lift at a given geometric incidence (Fig. 1). As the aspect ratio decreases, the 
upwash at the centre from the trailing vortices is gradually reduced, and there is a 
downwash over the whole wing span if the aspect ratio is small enough. None of 
these changes in the pattern of the streamwise vortices affect the chordwise loading, 
apart from an overall factor, as can be seen from Fig. 3. The same behaviour has 
been found on unswept wings.* The essential effect of sweep is manifested as a 
change in the pattern of the spanwise vortices and hence of the chordwise loading, 
rather than in any difference in the behaviour of the trailing vortices. 


The present theory is based on the supposition that the main effect of sweep 
can be deduced by considering the changes occurring in the downwash from the 
spanwise vortices, whereas the downwash from the trailing vortices is calculated 
as for a straight wing, i.e. from.a consideration of the vortex distribution in the 
wake. The effect of a kink in the spanwise vortices can best be studied by 
considering a swept wing of infinite aspect ratio with uniform load distribution, 


i.e. the case where there are only spanwise vortices of constant strength along the 
span. This is done in Section 2, which is a continuation of the investigations 
described in Ref. 1. The “second problem” of the classical aerofoil theory, that 
is to find the spanwise load distribution for a wing of given plan form, is treated 
in Section 3. Wings with given load distribution (“first problem ”) are dealt with 
in Section 4, and wings with minimum induced drag (“third problem”) in Section 5. 


Notation 


x,y,z rectangular co-ordinates measured in terms of the wing chord; 
x=0 at leading edge; x-axis in wind direction; z-axis downwards 


wing chord 

wing span 

spanwise co-ordinate, in terms of the semi-span 
aspect ratio 

chordwise position of the aerodynamic centre 
geometric angle of incidence 


effective angle of incidence, i.e. downwash induced by spanwise 
vortices 

induced angle of incidence, i.e. downwash induced by streamwise 
vortices 

angle of sweep of the mid-chord line 


*Effects introduced by small aspect ratios are not considered here. 
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FALKNER 


Fig. 1. Experimental and theoretical spanwise lift distributions and positions 
of the aerodynamic centre. 


velocity of the main stream 
induced velocity in z-direction 
downwash induced by spanwise vortices 
downwash induced by streamwise vortices 
strength of spanwise vortices along chord 
non-dimensional circulation along span =C,c/(2b) 
constant in equation (3) 
pressure coefficient [=(p— p,)/GpV,”)] 
difference of the pressure coefficients on upper and lower surfaces 
of the wing at any point 
local lift coefficient of section 
overall lift coefficient of the whole wing 
pitching-moment coefficient 
two-dimensional lift slope (usually assumed to be 27) 
sectional lift slope (=C,/,) 
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(1) VORTICITY VECTOR 
(2) SPANWISE VORTEX aot 
(3) STREAMWISE VORTEX 


Fig. 2. Vortex pattern near the centre of a swept wing (schematic). 


n 


an index, dependent on ¢ and y, defining the chordwise distribution 
of vorticity (see equation (3)) 


static pressure 

free stream static pressure 

air density 

aerofoil thickness 

interpolation function, used in equation (18) 


parameter in equation (26) which takes account of aerofoil 
thickness 


parameter in equation (26) which takes account of the lift reduction 
due to the boundary layer 


parameter used in equation (i) 
co-ordinate of aerofoil centre line at centre section of wing 
co-ordinate of aerofoil centre line at sheared section of wing 


2. The Chordwise Loading 

Consider the centre section of a swept wing of infinite aspect ratio, which is so 
shaped as to give a uniform loading along the span. The downwash v,, induced 
by the vortex system which replaces the wing can be written as: 
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Fig. 3. Experimental and theoretical chordwise loadings (the calculation includes thickness 
effects on the basis of an unpublished method by J. Weber). 


0-5 


as suggested in Ref. 1. The parameter o in this relation is given by 


1 
Yx (x) 
in the notation* of Ref. 1. Numerical values of « have been evaluated for various 
aerofoil sections of different thickness/chord ratios and several functions y, (x), 
giving the important result that « depends primarily on the angle of sweep ¢ and 
to a lesser extent on x and z. Ifo is written as a series of functions beginning with 
the term 


c= 


I, (x, z)— (¢) 


o=T tang, 


then the higher-terms, which depend on the section shape z(x) and ¢, are small 
compared with the first term. Thus, equation (1) can be replaced by 


1 
} 


*It may be noted that a system of co-ordinates is now used which differs from that in Ref. 1. 
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to a first order. This relation can be used either to find the downwash and hence 
the shape of the section for a given load distribution (see Section 4), or to determine 
yz(x) for a given aerofoil shape and incidence. In the latter case, it can also be 
assumed that equation (2) still represents a good approximation for the downwash 
at the centre section, even if the chordwise loading at the centre section differs from 
that farther out on the span. A typical example of this is the plane wing of infinite 
aspect ratio where v,/V,—,=constant everywhere, whereas y, varies with y as 
well as x. 


A solution of the integral equation (2), for v,/V,=2.=constant, which gives 
Yr—> © at the leading edge (x=0) and y,.=0 at the trailing edge (x=1) is, as 
suggested in Ref. 1, 


where C is a constant to be determined later and O<n<1. Inserting this into 
equation (2) gives 


{k, (n)-+[K, (n)+x tan 9] } 


by equation (8) of Ref. 1, where 


K, (n)= 


ds 
nJjs/"+1 sintn 
0 


co 


1 dt 
K,(n)= ~~ tan mn 


It follows that for symmetrical aerofoils, where va/V,=2, and so is 
independent of x, 


T 
- = =f tan 9. 
K, (n) tan =n ? 


This leads‘to a relation between the parameters n and 9:— 


n=3(1- 5), 


which was derived empirically in Ref. 1 but now follows directly from the down- 
wash equation (2). Further, 


C=2V,2, sin ™n 
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so that the chordwise vortex distribution for a plane wing becomes 
1-x\" 


or ,2, 008 »(4—*) 
since sin =2=cos¢ by equation (5). 


The chordwise load distribution can be expressed by the difference, AC, (x), 
between the pressure coefficient on the upper and lower surfaces of the aerofoil. 


From the Kutta-Joukowski theorem 


4 pV," Vo 


at the centre section where the direction of the vortices is normal to the free 
stream. Hence, 


AC, (x)= — 42, sin = 4,cos ¢(4—*)’ (8) 


The sectional lift coefficient, C,, is obtained by integration : — 


1 1 
| AC, (x) dx= 42, cos ¢ | dx, 
0 


0 


1 


and, since dx=nK, (n), 
0 


Cos ¢ 
this gives —— 
sin =n 


which may be expressed as 
Cy =4nrn 


or (1- a. 
The sectional pitching moment coefficient C,, about the leading edge is 
1 1 
AC, (x) x dx= — 42, sin =n x dx 
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which gives 
Cy= — [1 (=) 


1 
since xdx= 5 (n-1)K,(n-1), 


The position of the aerodynamic centre x,,, is then 


With this, the chordwise load distribution at the centre section and its properties 
are fully determined. 


At the sheared part of the wing, where the flow can be assumed to be two- 
dimensional, the downwash is again given by an equation of the type of equation 
(1), but with «=0, and its solution is of the type of equation (3). Hence, 


which gives K,(n)=0; ie. n=}, 


Be 


for a symmetrical aerofoil, and C=2V,z.._ The chordwise vortex distribution is then 


and the chordwise load distribution is 


AC, (x)= —2c08 
Vo 
instead of equation (7), since the vorticity vector is inclined at an angle ¢ to the 
main flow and only the velocity component V,cos¢ produces a lift force on the 
vortices. Thus 


—x\* —x\t 

AC, (x)= —42,c08 (4=*) = 42, c08 ¢ (4+=*) » 83) 
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To determine the chordwise loading at any intermediate station between the 
centre and the sheared part, use can be made of the fact that the expression for 
AC, (x) at the centre, equation (8), has the same form as that at the sheared part, 
equation (13). Only the values of the exponent n differ. It is therefore assumed 
that 


1-x\" sin 1-x\" 
AC, (x)= — 42, cos ¢ (+24) ) ‘ . (16) 
everywhere, which gives 
COs 9 
C.=4tn es : (17) 


as in the two limiting cases, equations (9) and (14). 
To obtain the variation with y of the parameter n, the following supposition 
is made: — 


This includes the case y=0, equation (5), if A=1 there, and the case of the sheared 
wing, y=00, if A=0 there. The interpolation function A (y) is related to the shift 
of the position of the aerodynamic centre 


A t= 
by equations (10) or (15). With equation (18), 


Since this is a physical problem, the locus of the aerodynamic centre along the wing 
span can be expected to be a continuous line which crosses the centre line at right 
angles; a reasonable approximation to this locus is the hyperbola whose asymptotes 
are the quarter-chord lines of the two wing halves and which passes through the 
aerodynamic centre at the centre section. This gives 


y being measured in terms of the wing chord. Since the term tan ¢/¢ does not vary 
very much with sweep, this may be simplified into 


A(y)= ov {1+ (2xy)?} ‘ 


These relations will be sufficiently accurate in most practical cases since it follows 
from equation (19) that an error < in the value of A causes an error of only </8 in 
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the position of the aerodynamic centre in the usual sweep range, around ¢=45°; 
the maximum error in x. can be ¢/4, if -=90°. From the relations (16), (17) and 
(19), the chordwise load distribution and its properties at any intermediate station 
between the centre and the sheared wing are fully determined if the parameter 7 
is taken from equation (18) with the interpolation function from equation (21). 


The tip of a swept wing can be regarded as part of the centre of a wing of 
opposite sweep. A consideration of the spanwise vortices leads to a downwash 
equation of the same type as equation (2), with a solution similar to equation (6), 
so that equations (12) and (18) are again valid, A being negative in the tip region. 
A more detailed investigation of the flow about the tip is not usually needed, since 
the lift falls to zero there. 


3. The Spanwise Loading 


The strength and distribution of the vortex system which replaces the wing 
may be determined from the condition that the downwash v, induced by the system 
is equal and opposite to the normal velocity component, «V,, of the main stream 
at the wing surface, thereby making the wing into a stream surface. Hence 


(x, y) 
Vo 

for a symmetrical aerofoil. The downwash v, is determined in two parts, 
Vz=Un+Vxz. Un is the downwash induced by the spanwise vortices which produce 
the lift force on the wing; it has been determined in the preceding section in the 
form v,=2V,. vz is the downwash induced by the streamwise vortices on the 
wing and in the wake, which do not contribute to the lift; its mean value can be 
written in the form v.=%V,. Thus the boundary condition reads 


The present theory is based on the suppositions that the aspect ratio of the 
wing is large and that sweep has no appreciable effect on the relation between , 
and the load distribution. In this case a is given by 

+1 
dy 


as in Prandtl’s classical aerofoil theory’ where 


(24) 


b/2 


is the non-dimensional spanwise co-ordinate and 


is a non-dimensional loading coefficient. 
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The value of z, on a swept wing of finite aspect ratio is also a function of the 
local lift coefficient, 


(25) 


by equation (17), where a is the “ sectional lift slope.” 


The relation gives the two-dimensional lift slope a, as 2 (for ,=0; n=1/2), 
which is true only for the two-dimensional fiat plate in inviscid flow. For two- 
dimensional aerofoils of non-zero thickness in viscous flow, the approximation 


a,=k(1+¢)2" ‘ ‘ (26) 


can be used, where « takes account of the thickness of the aerofoil and to a 
sufficiently close approximation «=0-8t/c. For swept wings, <=(0°8t/c)/cose. A 
measure of the lift reduction due to the boundary layer is given by k; k=1 in 
inviscid flow, k~0-9 for sections of about 10 per cent. thickness/chord ratio at 
Reynolds numbers around 10°, at small lift. A generalised relation for the sectional 
lift slope is thus 


a=a, COS ¢ (27) 


sin 

It remains to be decided how this relation for the sectional lift slope can be 
applied to wings of any given plan form. A simple way is to work on the hypo- 
thesis that in determining a it is essential at each spanwise station to represent 
correctly the mean sweep of the spanwise vortices and their chordwise distribution, 
and that it is relatively less important to represent correctly the conditions at the 
neighbouring stations. Thus, the conditions at the centre section of a given wing 
can be taken from a hypothetical wing (A) in Fig. 4, which has the sweep ¢.,2 of the 
mid-chord line of the given wing and infinite aspect ratio; the conditions at the 
sheared part of the wing can be taken from a hypothetical sheared wing (B) of 
infinite aspect ratio; and the conditions at the tip from a hypothetical wing (C) 
of semi-infinite span. The chord of each of these wings is equal to the local chord 
at the station considered. 


Equation (27) can be used for the centre and for the sheared part of the wing 
and also for any intermediate station. The hypothetical wing has its central kink 
at the same position as the real wing and the distance y of this station from the 
centre must then be measured in terms of the local chord at the station when 
determining the value of A(y) from equation (20). Equation (27) may also be 
applied to the tip region but in this case A, which is needed to determine n from 
equation (18), must have negative values, i.e. A= —1 at the tip. The value of A 
can be found from equations (20) or (21) by taking y as the distance of the station 
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Fig. 4. Replacement of a tapered wing of finite aspect ratio by hypothetical wings. 


considered from the tip, measured in terms of the local chord, and by reversing the 
sign. In such cases where the centre and the tip effects overiap, the sum of the 
two values of A may be taken. 


The boundary condition (23) can now be written 


+1 
1 (a) ar 
by equations (24) and (25), 
2 
dr’ 


a(n)c 
1 

This is the integral equation of Prandtl, but regarding the sectional lift slope a as 

a function of 7 according to equations (27) and (18). The most suitable method 

for solving it is that of Multhopp”’. Having found the spanwise loading from this 

equation, the chordwise loading can be determined from equation (16). 


The calculation of the load over the surface of a wing of given plan form is a 
matter of one or two computer hours. Some examples are shown in Figs. 1 and 3. 
The results agree well with the experimental points and the present theory 
adequately explains the peculiarities of the load distribution over swept wings. 
Small deviations between calculation and measurement are just apparent for the 
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wing of aspect ratio 3. This, then, appears to be about the lowest value of the 
aspect ratio to which the theory in its present form can be applied in the usual 
sweep range around ¢=45°. 


Some results from the calculation methods of Weissinger and Falkner“: *) are 
also shown. Both methods tend to over-emphasise the centre effect. 


4. Wings with Given Load Distribution 


A wing of given plan form can possess a given load distribution over its surface 
at one angle of incidence if it is suitably cambered and twisted. Consider a swept 
wing of infinite aspect ratio, as treated in Section 2. Let 


AC, (x)= - . 


be the required chordwise load distribution and let C,, be constant along the span. 
There are only spanwise vortices of constant strength along the span (i.e. «=z<,), 
and their chordwise distribution is 


2 cos ¢ COS 9 


by equation (12), because sheared wing conditions are required everywhere, the 
centre included. 


At the centre, the downwash equation (2) must be satisfied, 


= 


+tane(*=*)'} 


With the streamline condition v,/V,—dz/dx, the co-ordinates z,(x) of the centre 
line of the aerofoil section at the centre of the wing can be obtained by integration : — 


(x)= 


if z.(0)=0. 


E —sin-* (1-2x)+ | }. (31) 


COS 9 


At the sheared part of the wing, the aerofoil is a flat plate, so that 
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where _ Cu 
cos 9 

is the angle of incidence there. The difference in the shape of the chord line between 

the two stations is thus: 


Az astan 9 —sin-* (1 —2x)+ {1-1 | 


This then represents the necessary modification of the centre line in order to obtain 
the flat-plate distribution at the centre of the wing. It can be divided into 
two parts: — 
T 
which represents a twist, and 


Az (x)—Aax =} 2,tan ], (35) 


which represents a camber line. In this expression for the camber line, 2, (and thus 
C,) and ¢ appear only as factors multiplied by a function of x so that the camber 
lines are similar to each other under all conditions. This camber line has its 
position of maximum camber at 29 per cent. of the chord. The camber at the 
centre of a swept-back wing is negative (z is positive downwards). The twist is in 
the sense of increasing the incidence at the centre. 


For a wing of finite aspect ratio, another angle of twist is needed to compensate 
for the different induced angles of incidence 2, from the streamwise vortices. The 
value of a, for given C,(n) and c(n) can be determined from equation (24). In the 
relations just given, 2 must then be interpreted as the effective incidence 2,. 


An untapered 45° swept-back wing of aspect ratio 5 with camber and twist in 
the central region according to equations (35) and (34) has been tested and the 
results are plotted in Fig. 5. It will be seen that the lift distribution at the centre 
can be restored to that of the sheared part of the wing by the foregoing method. 
At the same time, this test confirms the adequacy of the downwash equation (2) 
which is the basis of the calculation. 


Another method of obtaining a given spanwise load distribution is to modify 
the plan form of the wing. Let C.=C,=constant be the value of the lift coefficient 
wanted at a given incidence z. The value of C,/ must be smaller than a,=a, cos ¢, 
and the difference between C,/a and a, will determine the aspect ratio of the wing. 
The chord distribution c(n) can be obtained from equation (28):— 

a(n) 
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CENTRE SECTION OF WING WITH — I-2 % CAMBER 
AND 3-1° TWiST 


BASIC SYMMETRICAL SECTION, t/c =0-12 
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Fig. 5. Experimental chordwise loadings on a symmetrical wing and on a 
wing with camber and twist. 


This is of a similar type to Prandtl’s aerofoil equation and can be solved in the 
same way (e.g. by applying Multhopp’s method). Since the value of a(n) depends 
on y/c and thereby on the unknown chord, the solution can only be obtained by 
successive approximations. The chordwise load distribution of such a wing is, of 
course, not the same at all spanwise stations unless the wing is both cambered and 
twisted. An example of the type of plan form obtained is shown on the right hand 
side of Fig. 6 (B). The inner part of the wing shows a slight inverse taper. 


5. Wings with Minimum Induced Drag 

The induced drag of a wing is at its minimum for a given lift if the induced 
downwash 2, from the streamwise vortices is constant along the span. This holds 
for both straight and swept wings. Since the same relation for 2 is used in both 
cases, equation (24), the spanwise load distribution must always be elliptical, 


where C, is the overall lift coefficient and A the aspect ratio of the wing. The 
induced incidence is then 
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Fig. 6. Calculated planforms (A) for minimum induced drag, and 
(B) for constant lift. The dotted line is the elliptic chord distribution. 


(38) 
and the induced drag coefficient is 


(39) 


Since 2; is constant, c,=a—«, is also constant along the span. Hence C,=a(n) a. 
is not constant, but varies directly as the sectional lift slope a(n). Thus, 
C,=a, cos ¢ 2, on the sheared part of the wing, but it has a smaller value over the 
central region and a greater value near the tips. If the wing plan form is such that 
the sheared wing value can be used as mean value, then 


C. =a, CoS? &=a cos ¢ 
L 0 e 0 TA 


by equation (38), which gives 


a(n) 
a, 


(40) 


14. 20.008 ¢ 


TA 
to the first order. Equation (41) is the usual result for an unswept wing of large 


aspect ratio, but with the two-dimensional lift slope replaced by the sectional lift 
slope of the infinite sheared wing. 
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With swept wings the plan form to give elliptic span loading is not elliptical. 
The chord distribution c(n) can be determined from equation (37), where y is 
defined as 


(y= 


which can be written as 


Ls 


579 C. 579 a, COS 


by equation (40). Insertion into equation (37) gives 


Now 


is the chord distribution of an elliptic wing with the same aspect ratio A. Hence the 
departures from elliptic chords c, and constant local lift coefficient C;, are given by 


 a,cosg _ 1 (42) 


The type of plan form having minimum induced drag is shown on the left-hand 
side of Fig. 6 as a highly tapered wing (A). This agrees with the general observation 
that delta wings tend to give a nearly elliptical spanwise load distribution. 


These results are consistent with the condition of R. T. Jones‘ (1951) that, for 
the drag to be a minimum, the downwash in the “combined disturbance field” 
(obtained by a superposition of the disturbance field of the given wing and the field 
in reversed flow of a wing with the same plan form, cambered and twisted to have 
the same lift distribution in reversed flow) must be constant at all points of the 
wing surface. The plan form of equation (42) ensures that the downwash in the 
combined field from the infinite ribbon of trailing vortices is constant everywhere, 
at all values of C,. To satisfy Jones’s condition, however, it would appear that the 
minimum-drag wing should also be twisted and cambered so as to give a constant 
downwash from the combined field of the bound vortices. By equation (2), the 
centre terms +7 tan ¢y,(x) cancel one another. For the remaining term in equation 
(2) to give constant values over the wing surface, y,(x) must be the two-dimensional 
flat-plate distribution, equation (6) with n=4. This can be achieved only at one 
value of C,, and a modification of the wing shape as explained in Section 4 is needed, 
with twist and camber as given by equations (34) and (35). Thus Jones’s condition 
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leads to wings with elliptic spanwise loading and “ flat-plate” chordwise loading. 
The flat unswept wing of large aspect ratio and elliptic plan form naturally fulfils 
these conditions’. Swept wings require a different plan form and must be cambered 
and twisted. 
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Note on the Stabilisation of Long Struts by 
an Elastic Medium 


W. J. GOODEY, A.F.R.Ae.S. 


1. Introduction 


This problem arose in connection with the stabilisation of the spar flanges of 
an aeroplane, where the spar passes through the fuselage. The fore-and-aft 
stabilisation of a compression flange in these circumstances may often be effected by 
means of intermediate supports, which are necessarily flexible, and the reliability of 
any calculations of buckling load is somewhat doubtful, mainly due to the difficulty 
of estimating the stiffness of the supports. In view of this it was decided in a recent 
design to consider a continuous support for the spar flange in the form of a thin 
horizontal web riveted to the flange and to the sides of the fuselage, and having a 
small member attached to its other edge. This arrangement gives an interesting 
case of a strut on a continuous elastic support, in which the support is not of the 
simple kind usually considered, where the reaction is proportional to the deflection, 
but is of a more complicated type involving interconnection between the supports. 


Notation 
Ox, Oy axes parallei and perpendicular to the strut, the origin being taken 
at the mid-point of this member in its undeflected position (see 
Fig. 1) 
2/=distance between points of application of end load to the strut, - 
which will be assumed to be simply supported at these points 
width of stabilising sheet, referred to as the web 
thickness of the web 
k,t effective thickness of web+reinforcing members, for taking end 
load in the y-direction 
N.,Ny,Nzy — stress loadings (i.e. stress x effective thickness) in the web and 
reinforcing, the first two being tension parallel to Ox and Oy 
respectively and the last being shear. The sign convention for shear 
is shown in Fig. 1 (N,=0 by the assumption made later) 
u, v components of displacement of the point (x, y), parallel to Ox and 
Oy, respectively 
Ua, Va 
upg, Up 
M., 1., Aa bending moment, moment of inertia about its own neutral axis, 
and cross-sectional area, of the strut (member AB in Fig. 1, alter- 
natively referred to as member 2) 


First received March 1952. 
[The Aeronautical Quarterly, Volume IV, August 1953] 
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corresponding quantities for the edge member CD (or member /) 


compression in AB 
tension in CD 
“Euler” buckling load for the complete beam ABCD considered 
as a strut simply supported at each end=*EA.A gh? /[L? (A.+ Ag)] 
Young’s modulus 
modulus of rigidity = E/[2 (1+.)] (assume strut and stabiliser to be 
made of the same material) 
Poisson’s ratio 

F,G,HandK _ functions of x, introduced in Section 2 


(P.=P, Ps=0 at x= +1) 


D_ denotes the operator 


A,(r=1,2,3) the roots of the cubic equation (14) 
n any integer 
m_ any odd integer 
w deflection of the point (x, y) of the web, normal to the plane Oxy 
Crm Cyy Cry tensile strains in x- and y-directions, and shear strain, respectively, 
in the middle plane of the web (i.e. mean strains in the web) 


2. The Distribution of Stress in the Sheet 


The stabilising sheet will be assumed to be reinforced with stiffeners perpen- 
dicular to the strut. The stiffener-sheet combination will then be capable of taking 
shear loading N,,, and end load in the y-direction N,, but its capacity for taking end 
load in the x-direction N, will be negligible for practical purposes, except in the 
immediate vicinity of the edge members. These assumptions lead to a comparatively 
simple solution, and are justified by the consideration that the buckling strength of 
the web in practice is likely to be negligible compared with the strength of the strut. 


The simplest approach would clearly be to consider the beam ABCD as a strut 
with an offset end load, using the known formulae applicable to this case. It was 
chiefly in order to find out how much in error the simple theory was likely to be 
that this investigation was originally made. The equations of equilibrium of a 
rectangular element of the stiffened sheet, assuming N, to be zero, are 

_ 


dy =0 and 


ON, _ 


+ =0 ‘ (1) 


The first of these shows that N,, is a function of x only, or 
Substituting this in the second of equations (1), and integrating with respect to y 


then gives 
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P= APPLIED 
END LOAD 


WEB STIFFENERS 

DIRECTION 
STABILISING WEB 


MEMBER 
Fig. 1. 


dF 
Ny=G-y7 


G being a function of x only. 


3. The Stress-Strain Equations for the Sheet 


On the assumption that considerable buckling of the sheet takes | place, 
expressions for the strains e,,., ey, and e,, Must be used which take into account the 
deflection of the sheet normal to its original plane. 


These are 
Ou 


ow dv dw \? 
Ox +4 » w= +4 


and _ ou dw Ow 
Oy ax + Ox dy 


(4) 


It is now assumed that the strut AB is completely stabilised in a plane at right 
angles to the sheet, and that the web stiffeners and the member CD remain straight. 
Buckling will then take place approximately in the x-direction, which means that w 
will change only slowly with respect to y, except near the edge members AB and 
CD, and even here, the derivative dw/dy is likely to be small enough for its square 
to be negligible compared with dv/dy. On the other hand, the maximum value of 
dw/dx may be quite large, particularly if the web stiffeners are closely spaced, but 
the product (dw/0x)(dw/dy) will remain practically zero everywhere, since when 
Ow/dx is large dw/dy will be negligible, while when dw/dy is appreciable (near AB 
or CD) dw/ dx will be negligible. 
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Thus, the only term in w which need be retained in equations (4) is the term 
4 (0w/ 0x)? in e,.. 


The other two strains are thus expressible in terms of derivatives of u and v 
only, and the associated stress-strain equations are then 


ov dF 


Ou 


and Ct 


dv 
+5) =Ny=Fo © 


Equation (5) may be integrated with respect to y, giving 


Ek,tv=H + Gy -4y° 


where H is an arbitrary function of x. 
Substituting this result in (6) and again integrating with respect to y, gives 
d°F 


dH 4G, 


where K is another arbitrary function of x. 


Equations (2), (3), (7) and (8) constitute a general solution for the mean stresses 
and the displacements in the plane of the stiffened sheet, on the given assumptions. 


4. The Boundary Conditions 


The particular form of the functions F, G, H and K depends on the conditions 
obtaining along the edges of the sheet. The equation for the deflection of the 
member AB is . 

d‘v, 
dx* 


d’v, 


El. dx* 


+P, +N,=0, 


this being the general equation for a uniform strut. 
Now the end load in AB is EA, 6u./ 0x (tension), and the shear loading F in the 
adjacent skin is equal to the rate of increase of this end load with x. 


Hence F=EA, 


The values of v., u. and N, are obtained from the general solution by taking y=0. 
Substitution in these two equations then gives 
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and 


where D denotes the operator d/ dx. 


The equilibrium of the member CD, for which y= -h, is now considered. 
It is assumed that it has practically no bending stiffness (/;=0) and an end load 
which is small compared with P,. The lateral loading (N,) along this edge of the 
sheet will therefore be negligible, and hence, from (3), 


The end load in CD is EA Oug/0x (tension), and the rate of increase of this quantity 
with x must be equated to - F, 


= =0 


ie. F+EAs; 3x 


Substituting for us from (8), in which y= - h is taken, gives 


F+ [K - {2k,(1+o0)F- DH}h-4h°DG-23h°D*F]=0 (12) 
y 
Equations (9), (10), (11) and (12) are four simultaneous differential equations 
with constant coefficients, provided that P. in (9) is taken to be constant. In fact, 
P, is equal to a constant plus a small variable part; this is obvious in considering 
the deflection of the beam ABCD from its unloaded position. 


The elimination of K and G from (9) and (12), using (10) and (11), then gives 
two relations between H and F. Finally, the elimination of F between these two 
leads to a differential equation for H, and therefore for v.. Omitting the inter- 
mediate steps, which are quite straightforward, the result is :— 


{ ( +. + 2k, (1 +0) AD? + { +P.D*} 4 | H=0 
(13) 


Considered as an equation for D*H, this has a general solution of the form 
(a, sin A,x +b, cos A,x), r= 2,3. 


where A,”, A,? and A,” are the roots of the cubic in A? 
{ kat (4 + +2ky (1 +0) ha? + { P.- ELM} —Ek,th=0 
D’H is now integrated twice, with respect to x, bearing in mind that H =Ek,tv., 


and that v.=0 at x= +1, hence obtaining the result 
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i= Ek,tv.=%,{ (+ sin sin + cos 2.x) } 


5. The Buckling Load for Pin-Jointed Ends 


The critical buckling condition for the case where the member AB is pin-jointed 
at each end is now investigated. This appears to be the only case which leads to a 
simple expression for the buckling load, and in fact for this purpose it is necessary 
to make the additional assumption that the points A, D and B,C are tied together 
by inextensible members AD and BC, so that the displacement v is constrained to 
remain zero at x=+/. In practice, there will probably be a certain amount of 
fixation of AB at each end, while the condition v=0 can only be satisfied approxi- 
mately, but as these departures from the assumed conditions will tend to have 
opposite effects on the buckling load the result of the present investigation sould 
be a useful guide for design purposes. 


The condition of pin-jointed ends means that 0°v,/0x*, and hence D?H, must 
be zero at x= +/, and hence from equation (15) 


>,a, sin cos A,J=0 . (16) 


The displacement v may be expressed in terms of H and its derivatives by the 
elimination of G and DF between equations (7), (9) and (11), whence 


=Ek,tH y (1+2 (P.D*H + 


At x= +l, v=H=D°H=0, and so D‘H=0. 
Therefore sin A,l=%,b,A,? cos A, =0 


The remaining requirement is that P, =P, Ps=0 at x= +1. 


Now 
Oup 
and Ps= + EAs ax 


while differentiating equation (8) with respect to x, and making use of equation (11) 
to get rid of F, it is found that 


Ekyt ~ { 2k, } -y (4+ 2 


Taking y=0 in this equation gives 


DK =Ekyt = — | 
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and putting y= — A in the same equation gives the result 


h* 
kt (Zt + 42) = 3 DG, 


G and D’G being given by equation (9). 


Therefore 


and from equation (15) this easily leads to 


sindA,l= 0 
(18) 
Taking the first of each of the pairs of equations (16), (17) and (18) it is seen that 
a, =a,=a, =0, unless one of the A’s is equal to nz/I, in which case the corresponding 
a, is arbitrary. 


The constants b, may be found by solving the second of each pair of equations, 
unless one of the A’s is equal to mz/(2/), where m is an odd integer, in which case 
the corresponding 5, is arbitrary. 


Thus, if A=nz// is a root of equation (14) the deflection of ABCD will be made 
up of a symmetrical component of a definite amount together with an arbitrary 
unsymmetrical component, while if A= mz/(2/) is a root of equation (14) the deflec- 
tion will be wholly symmetrical but of arbitrary amount. Both conditions are 
covered by the equation A=nz/(2/), where n is any integer; i.e. A=nt/L, where 
L=2l. Putting this value of A into equation (14) gives 


(19) 


Buckling will occur with the integral value of n which gives the lowest value of P 
from the formula. In practice the stabilising sheet will probably be designed so 
that this happens when n=1. 


The first term in equation (19) is the buckling load for the member AB 
unsupported, except at its ends. Taking only the first term inside the square 
brackets gives the “Euler” buckling load P. for the beam ABCD considered as a 
strut simply supported at each end, and neglecting shear deflection of the web and 
flexibility of the lateral stiffeners. The other two terms inside the brackets represent 
the effect of the web and stiffener flexibilities. 
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TABLE I 


(tons) 


312°9 
335°5 
353°4 
367:0 
383°3 
388°5 
385-9 
378°5 


It will usually be found that the first and last items in equation (19) can be 
neglected in comparison with the other two. The buckling load is then given with 
n=1, leading to the simplified formula 


P=P.(1+ 


which is essentially the result given by Timoshenko in his Theory of Elastic Stability 
(Art. 26). 


6. Numerical Example 
Taking A,=2:0 in.*, As=0-2 in.*, t=0°036 in., L=SO in., ky=4, 
E=4,690 tons/in.?, e=4, and leaving the depth A variable, equation (19) reduces to 
0:285123 


P=6-172n? + 18-5154 F 


+0-00043293 hn’ | tons, 


with f in inches. The value of k, has been assessed on the basis of 4 in. x 4 in. ~ 
0-036 in. stiffeners on the web, at 6 in. pitch, taking 1 ‘n. of web with each stiffener. 
With this pitch of stiffeners the buckling load calculated for a continuous distribution 
of stiffening material should be reasonably accurate. 


The value of P has been calculated for various values of h from 8 in. to 20 in., 
and with n=1, 2, 3, 4, 5 and 6. The results are given in Table I. 


The lowest critical loads for each value of h are in italics. The best value 
obtainable in this case is 270 tons at a web depth of about 11-6 in., and is about 
60 per cent. of P.. 


in) n=1 | 2 | 3 4 | 5 | 6 | (tons) 
8 | 1543 | 2856 | 298-0 | 2908 | 301-7 | 3360 | 2155 
9 185-9 | | 304-8 | 287-2 | 204-8 | 328-7 | 2727 
10 «(217-0 | 307-2 | 281-7 | 2876 | 321-8 | 3366 
11 249-4 | | 306-3 | 275:2 | 2804 | 3154 | 4073 
12 | 281-9 | | 303-2 | 2682 | 2736 | 3096 | 48488 | 
14 | 3460 | | 2928 | 254-2 | 2613. | 2996 | 65988 
16 | «406-7 | | 279-8 | 241-2 | 2508 | 2913 | 
18 =| 4621 | | | 2295 | 2419 | 2846 | 109079 
20 | | | 253-1 | | 2344 | 2789 | 134668 ™ 
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Supersonic Flow Past Wing-Body 
Combinations 


W. CHESTER 


(Department of Mathematics, University of Bristol) 


Summary: The supersonic flow past a combination of a thin wing and a slender 
body of revolution is discussed by means of the linearised equation of motion. The 
exact equation is first established so that the linearised solution can be fed back and the 
order of the error terms calculated. The theory holds under quite general conditions 
which should be realised in practice. 


The wing-body combination considered consists of a wing symmetrically situated on 
a pointed body of revolution and satisfying the following fairly general conditions. The 
wing leading edge is supersonic at the root, and the body is approximately cylindrical 
downstream of the leading edge. The body radius is of an order larger than the wing 
thickness, but is small compared with the chord or span of the wing. 


It is found that if the wing and body are at the same incidence, and the aspect 
ratio of the wing is greater than 2 (M®— 1)~!, where M is the main stream Mach number, 
the lift is equivalent to that of the complete wing when isolated. If the wing only is at 
incidence then the lift is equivalent to that of the part of the wing lying outside the body. 


The presence of the body has a more significant effect on the drag. If, for example, 
the body is an infinite cylinder of radius a, and the wing is rectangular with aspect 
ratio greater than 2(M?-—1)-!, then the drag of the wing is decreased by a factor 
(1—2a/b), where 24 is the span of the wing. 

When these conditions do not hold the results are not quite so simple but are by 
no means complicated. 


|. Introduction 


The linearised theory for supersonic flow past thin wings and slender pointed 
bodies is now well established and has considerable literature (see References, 
particularly the extensive bibliographies in Refs. 6 and 13). The aerodynamic 
forces predicted by this theory agree with experiment sufficiently for the theory to 
be used by aerodynamicists as a basis for calculations. Recent experiments, how- 
ever, have shown that the drag of a combination of wing and body differs by as 
much as 20 per cent. from the sum of the drags of the wing and body when 
separated. In practice linear theory is more accurate than this, which suggested 
that at least a partial explanation of this interference drag may well be obtained 
from a theory based on the linearised equation of motion. 


The possibility was further impressed on the author by the work of Lord"? 
on flow past slender bodies, where use was made of the fact that the perturbation 
velocities in the vicinity of a slender body, whose radius is small and O(a), can 


First received December 1951. 
[The Aeronautical Quarterly, Volume IV, August 1953] 
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be obtained on a linear theory with an error which is at most O(a loga). Thus, 
boundary conditions on the body which are larger than this, say O(a’), are 
significant. The author would like to acknowledge the fact that it was the signi- 
ficance of this result which led to the present investigation. 


The wing-body combination considered consists of a pointed body of revolu- 
tion having a nose, defined as that part of the body upstream of the wing leading 
edge, for which the cross-sectional radius R and its derivatives in the direction of 
the main stream are small and O(a) (the unit of length is taken to be the wing 
chord at mid-span). Downstream of the leading edge, R=a+O(a’) and its 
derivatives are O(a’), so that this portion of the body is to the first order a cylinder 
of radius a. The wing ordinate is uniformly small and O(n), where » is the maximum 
thickness, the wing leading edge is supersonic at the root and the span is at leas: 
comparable with the chord at mid-span. That is to say the results are not significan 
for slender wing-body combinations; there is, indeed, an alternative theory for such 
bodies"* *7), 


The formula obtained for the drag consists of the contribution arising from the 
potential fields of the wing and body separately, which is O(y?+a*) at zero 
incidence, together with a further term whose magnitude is O(an?). The error 
is O (n° + na‘ log” a+ log? a). If a and are comparable then, of course, the 
extra term is no more significant than the error term. In practice, however, a is 


larger than 7; in fact a survey of one or two typical wing body combinations indicated 
that »=O (a’) is nearer the truth, in which case the interference effect is significant. 


There is no interference effect on the lift at incidence if the wing is symmetrical 
and at zero incidence relative to the body axis; within the approximation the lift is 
effectively that of the wing itself. If the wing is cambered, then there is an extra 
term which is O(an) compared with the wing contribution which is O (y+) if the 
angle of incidence is ¥. 


Notation 


velocity of sound 

approximate radius of the body downstream of the leading edge 
semi-span of wing when rectangular 

pressure coefficient 

distance from the nose of the body to the leading edge of the wing 
component force in the direction of the body axis 

see equation (47) 

wing ordinate relative to maximum wing thickness 

Bessel functions of the second kind 

component force perpendicular to the body axis 

length of body 

Mach number 


aN 


<a 
a 
NY 
288 


SUPERSONIC WING-BODY COMBINATIONS 


in equations (4), (64) and (66) p denotes the pressure, otherwise it is 
the Heaviside operator 

see equation (3. 

local speed 

radius of body 

polar co-ordinate 

cross-sectional area of body 

variable of integration except in Section 7, where it is an element of 
length of the curve 

variables of integration 

main stream velocity 

Cartesian co-ordinates 

(M? — 1) 

Euler’s constant in equations (18) and (43); otherwise adiabatic 


index 
maximum thickness of the wing 


polar co-ordinate 
density 
spanwise distance from tip of wing (see Fig. 2) 
velocity potential 
velocity potential of perturbed field for body alone at zero incidence 
additional term in the velocity potential required when the wing is 
added, the body being at zero incidence 
additional term in the velocity potential arising from incidence 

= Und,w+Uv4.~, velocity potential of the perturbed field of the wing 
alone 
interference potentials, the contributions to 9, in addition to 9,» 
interference potential, the contribution to ¢, in addition to 9. and 
the contribution from the body alone at incidence 
angle of incidence 
variable of integration 


Suffixes, etc. 


Suffixes “w” and “b” denote wing and body. Other letter suffixes denote 


differentiation. 
’ denotes an interference potential. 


+ and ~ denote the upper and lower surfaces of the wing. 


2. Formulation of the Problem 


Consider a slender pointed body of revolution equipped with thin wings having 
a chord at mid-span of unit length. Tie body is at rest in a uniform supersoni¢ 
stream of velocity U inclined at an angle ¥ to the axis of the body, 
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Upstream of the leading edge, the radius R of the body, and its derivatives along 
the axis, are O(a). Downstream of the leading edge R=a+O (a’) and its derivatives 
are O(a’). The maximum thickness of the aerofoil is ». 


Let O (x, y, z) be a right-handed system of co-ordinates, with origin at the nose of 
the body, such that the z-axis is along the axis of the body and the x-axis in the 
direction of the span of the wing. The wing is symmetrical about the plane x=0 
and lies approximately in the plane y=0. 


Polar co-ordinates r, 6, are defined by x=rcos@, y=rsin@, and the following 
expansion, in terms of the parameters » and y, is assumed for the velocity potential 
of the-resulting flow: — 


@=Uzcos¥+Uy sin U9, (z.r)+ Ung, (z, 7,6) +U 9) + O 


where 9, is the velocity potential of the perturbed field arising from the body alone 
at zero incidence, #, is the extra term resulting from the presence of the wing, and 
v(rsin 6+ 9,) is the extra term arising from the axis of the body being inclined to 
the direction of the main stream. : 


The velocity is then given by (cf. Ref. 11) 


+9." 


2 


= (1+ + (or + + ¥ Sin + + O(n? + 


QO, =2 (1 Po2) Piz 24 orP ir 


Q, =2 (1 Po:) Go: + (sin 6+ or). 
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When the motion is irrotational, the equations of momentum and continuity 
become . 


1 
V q*)+ —Vp=0 
(4) 
and (pVo)= 0. 


Elimination of the pressure p from these equations gives 
2A°V*o= (Vq’). (V9), 


where A is the local velocity of sound and 


where A, is the velocity of sound in the undisturbed flow and y is the adiabatic 
index. Substituting expansions (1) and (2) for ¢ and q in equations (5) and (6), gives 


2 { —4(y- +10, + vo.) } { + + WV } 
+ {Qor + + VOor} {bor + +¥ (Sin 4+ +O (4? + 
where M=U/A, is the Mach number of the undisturbed flow. 


Equating coefficients, 


2 (y- 1)Q, = Q,.(1 Qo: Por (8) 


tig -20-DQ } Qiz(1 + (9) 


1 
2 4 (y- 1)Q, } Vo. (y- 1)Q.V 70, = + Qo (1 Py2)+ 
+ (SiN 9+ (10) 
Note that when these equations are linearised, ,,%,,%, all satisfy the equation 


The boundary conditions for ¢ are: — 


(A) On the body 
=R,(z) when r=R(z), 


Port ir + (Sin + 012+ 
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{ 1+ oz} 
sin 9+ ¢2,= Rive: 
when r=R (z). 


(B) On the wing 
If the upper surface of the wing is defined by 


then the direction cosines of the normal are proportional to (—h.*,1,—h.+), and 
the boundary condition of zero normal velocity is therefore 


nh,* oz nh.*¢.=0 on y= nht 
or — on y=0. 


Substitution for ¢ from equation (1) yields for x>0, 


Thus the boundary conditions on the part of the plane y=0 occupied by the wing 
are, for x>0, 


1 
Tr — 
For negative values of x, the boundary conditions follow from symmetry. 


3. Solution for % 


A first approximation to ¢, will satisfy equation (11) subject to the boundary 
condition (see equations (12)) 


(1 + 
=j(z) (say) 


on r=R. 
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In Heaviside notation, with p~'= | dz and 2?= M?—1, the solution of equation 


0 


(11) having symmetry about the z-axis and behaving correctly at infinity is 
¢,=CK, (apr), ° ‘ (17) 
which reduces to 


when r=O (a), where C is a constant determined by the boundary condition (16) and 
y now represents Euler’s constant. In fact substitution of (18) in (16) gives 


C=—-—{1+0(a’)} Rf. 
so that —{1+0O(a’)} RfK, (apr). 
Since 9, is O(a? log a~') when r=O(a), it follows from (16) that 


f={1+O(@ loga—')} R., 
so that, finally, 
— {1+O(@ log a~')} RR.K, (apr) 


—ar 


2 (z—s)? — 


{ 


0 
where S=R? is the cross-sectional area of the body. 
This result was originally obtained by von Karman and Moore“. 


Having obtained the solution of the linearised equation, the magnitude of the 
second order terms is now calculated from the exact equation for ¢, (equation (8) ). 
From (21) it is deduced that 


in a finite region containing the body. If this is substituted in the second order terms 
of equation (8), the result is 


(23) 


The correct solution for ¢, may thus be written, with the help of (21), 


z-—aer 


be 1 S,(s) ds 
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provided that r=O(1). In particular this will hold on the surface of the wing. 
Beyond the wing tips ¢, is uniformly O(a’) and thus the error term here is O(a‘). 


4. The Accuracy of the Linear Solutions for ¢., 9. 


Assuming that ¢,,9%, and their derivatives are O(1), and substituting the value 
of , given by equation (22) in equations (9) and (10), it is found that both ¢, and 9, 
satisfy the equation 
O(a’) 


24 —M?2o..— 
—M’9.: 


(25) 
if r=O(1). A particular integral is O(a’). Since a particular integral for large r is 


also O(a"), the solutions for ¢, and ¢, based on the linearised equation (11) can be 
said to be uniformly accurate to within a term of this magnitude. 


5. Solution for % Upstream of the Leading Edge 
In polar co-ordinates, equation (11) becomes 
1 1 
Porr + Tr Por + r = (26) 


If a solution for ¢. of the form 
is assumed then the equation for g is 


and the boundary condition (12) becomes 


onr=R. 


The appropriate solution of (28) is 
g= BK, (zpr) ‘ ‘ (0) 


which reduces to g= {1+O(a@ loga™')} 


when r=O(a). The constant B is determined by the boundary condition (29) which 
gives, for ¢,, the solution 


apR°K,(apr)sin@ 
_ {1+ O(@ loga~')} sind S. (s)(z—s)ds 


(33) 


0 
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If r=O (a), this becomes 


r 


The solution was first given in this form by Lighthill®®. 
From Section 4 it is deduced that the correct solution differs from (33) by a 
term which is O (a’). 


6. Solution Downstream of the Leading Edge 


In this region the velocity potential 4. of the wing is introduced; this is the 
solution of the linearised equation for flow past the wing alone which, for con- 
venience, is assumed to be continued in some suitable manner into that region of the 
plane y=0 actually occupied by the body. For present purposes it is sufficient that 
the continuation be such that the wing is symmetrical about the plane x=0, and 
that the wing ordinate is continuous and O(»), so that spanwise variations are 
O (an) inside the body. 


Writing dw = + 
then =h,*+ | 
_ | 
roo 


on the upper surface of the wing. 


Remembering that R=a+ O(a’) in this region, we have (say) 
where, by (15), 
_ ht + 
r )or+ h, Poz 


G7) 
= (1+ 0(@)} K, (apr) 


for positive values of x on the wing. Also 9,’ and ¢,’ are symmetrical about the 
plane x=0. 
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The terms on the right of equations (37) are significant only when r is small, and 

for the present purpose these relations may be simplified to 
« 
Ly» O@) 


r 16 r? 


a’ 
P20 {1+ O(a)} r 


on y=0, x>a. 


Beyond the wing tip the condition is that 9,’ and 9,’ should be continuous. 
Where the span of the wing is large compared with a, no significant discontinuity 
will be produced by assuming that equations (38) hold for all r on y=0. If the 
leading edge is swept back there could be a region near the leading edge where a 
discontinuity will appear beyond the tip of the wing, but this will produce a 
correction term of the same order as the terms in ¢,’ and ¢,’ which arise from the 
boundary conditions (38). Now it will appear later that these latter terms do not 
contribute significantly to the aerodynamic forces acting on the wing-body com- 
bination. On the assumption that the terms required to correct an insignificant 
contribution to the potential will themselves be insignificant, the boundary condi- 
tions given by (38) will not be modified. It is necessary, however, to assume that 
the leading edge is supersonic, otherwise the term in 4,’ which arises from satisfying 
the boundary conditions on the body violates the condition of continuity on the 
plane y=0. This point is discussed further in Section 7. 


Boundary conditions on the body are now violated by the extra terms intro- 
duced to take account of the wing potential. Thus, by (12) and (36), 9,’ and 4,’ must 
also satisfy the conditions 


(39) 


Ou = — Piwr +O (a’) 
= — Powr +O (a’) 


on r=a. 


It is now assumed that ¢, is analytic near the body, so that we may write 


Vi = 
or oy ),_,sin@+0@ 

=(nh,.* sin (a) 


where 7h, is the wing ordinate in the plane of symmetry of the wing, x=0. 
Equations (39) then become 


= sin6+O(a) 


$27’ =sin 6+ O (a) 
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The boundary conditions (38) and (41) are sufficient to determine ¢,’ and ¢,’ 
uniquely in the upper half-plane. The complete solution is somewhat complex, but 
considerable simplification is effected by using the fact that the boundary 
conditions are significant only in the region of the axis. It is known that in this 
region derivatives of the potential in the z-direction are not significant compared 
with derivatives in the other two directions, and the potential behaves like a two- 
dimensional solution of Laplace’s equation in the variables x and y (see Ref. 17 for 
further details). This suggests that the solution of the present problem will be 
indicated by the solution of Laplace’s equation which is symmetrical about the 
plane x=0 and, for x >0, reduces to 


C(r) on y=0 for +>a 


on r=a 
(say). Such a solution can be constructed by a distribution of sources on y=0 and 
on r=a, in fact 
1 
[tog (r? +A? + 2Ar cos 4)! + log (r? + A* 2Ar cos 6)! 2 log r| 


2 [tog (r? +A? + 2Ar cos 6) + log (r? +A? 2Ar cos 4) ‘| 
A 
0 


+ Jew [og {r? +a? 2ar cos (6 + log {r? + a” 2arcos (4+ |ay-+ 


is the required solution. But, when p is small 


K,e)= - {1400} } 43) 


Thus it follows, from (42) and (43), and the boundary conditions (38) and (41), that 


|| siny || Ko {ap (7? 2arcos (6 - 4K. (apr) | dy + O(a" loga~') 


(44) 
Kol ap(r? +A? + 2Ar-cos + Ky{ 2p(r? + A? — 2Ar cos — 2K (apr) + 


+ {ap (r? +A? + 2Arcos 6)*} da | |siny (2 


~2ar cos (6—¥)}#} (apr) dy + O log a-") 
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for these relations reduce to the form (42) when r is small. Since they also satisfy 
the linearised equation (11) and have the correct behaviour at infinity, all the 
conditions are satisfied within the approximation. 


When interpreted, 9,’ and ¢,’ are given by 


2H a(r* +a" ar cos _ 
a= =| [sin dy | — a (r? + a? —2ar cos )}3 


H {z-s-ar} 


| +O(@loga") . . (46) 


H(z)=1 z>0 
(47) 


=0 z<0 


is the Heaviside unit function, and c is the distance of the leading edge from the 
nose, measured along the axis, 


co 
cosh a (P+ + 2Ar cos H{z—c—a(r?+A?+2Ar cos + 


+ cosh Arcos H {z-c- a(r?+A?—2Arcos | + 
+ _ | ddcosh a(P? +)? 4 2arcos 6)! H {z-c- a(r? +A? + 2Arcos - 


-a 


a 
“| |siny | dycosh- { 


x H {z-c- +a? -2arcos(—¥)!} + O(a? loga-") for y<0. . (48) 


The field for y<0 can be deduced from (46) and (48) by replacing /».+ by 
—ho.~ , if the lower surface of the wing is defined by 
y= nh- (x, 2), ‘ : . (49) 


and from the fact that ¢,’ is anti-symmetrical about the plane y=0. 


7. Further Modification of ¢%,’ 
It may be noted, from (44) and (45), that when r is large compared with a, 


= 2" h,.+K, (apr) + O (a? log a-") 
(@’ 
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The contribution of the first term in ¢,’ is restricted to the region lying inside the 
Mach cone emanating from the junction of the wing leading edge and the axis of 
the body. If this cone does not run off the wing before the trailing edge then, 
upstream of the trailing edge, the condition that the potential should be continuous 


across the plane y=0 beyond the wing is not violated to a significant degree and the 
complete solution upstream of the trailing edge is obtained without further modifi- 
cation. Even if the aspect ratio of the wing is sufficiently small for the Mach cone 
to run off the wing at the tip, the potential will still be continuous if the wing is 
symmetrical about the plane y=0, for then ¢,’ is also symmetrical about this plane. 


Thus the only further modification of the solution occurs when the wing is 
not symmetrical about y=0 and the Mach cone meets the wing tip at some point O’ 
upstream of the trailing edge (see Fig. 2). The solution can be continued down- 
stream of O’ if a potential field can be found which is confined to the Mach cone 
with apex at O’ and which annuls the discontinuity in potential produced by 9,’, but 
leaves 0¢/dy unaltered on the wing. To illustrate the procedure in this case the 
wing tip is considered to be parallel to the axis of the body downstream of the 
point O’. The problem can then be solved with the aid of the Green function for a 
semi-infinite barrier. 


It is known from Green’s theorem that if ¢ and G are two solutions of the 
equation (V?—?p?)¢=0, ¢ and its derivatives being continuous inside and on a 
closed curve C, and G being continuous except at a point P, where G~—logp, p 
being the radial distance from P, then 
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Cc 
where n is the inward drawn normal to the curve C. 
If, in addition, 


G=0 onC, 


When the curve C consists of the two sides of the positive x-axis together with a 
large circle whose radius tends to infinity, the Green function satisfying all the 
above conditions has been derived by Gunn” in the form 


G(r, 9; r,,9.=F {r,r,.(0-9,)} -F {r.r.,(0+9)},  ‘©3) 


where (r,,9,) are the polar co-ordinates of the point P, and 


2(rr,)tcos 


(r? +r,? 2rr, cos 


8. Application of the Green Function 


A consideration of the expressions for 9,’ in the upper and lower half-planes 
shows that they will produce a discontinuity of magnitude 


on y=0 beyond the wing tip downstream of the point O’. A further term ¢,”, which 
is anti-symmetrical about y=0, will annul this discontinuity if it satisfies the follow- 
ing boundary conditions in the upper half-plane: — 


=~ hye )K, (ap y=0, >0 


20,” 


where b is the semi-span and the origin has been transferred to O’, so that 
x=xX +b. 


300 


| 
= 


1) 


A) 


at 


SUPERSONIC WING-BODY COMBINATIONS 


In terms of the Green function just described, it is found that 


9,” 6.)= = (hys* hy:~) x 


» {p(x +b)} 2 


(57) 


where F’(-96,))= (x’, ro’. P (58) 


ay 
and (r,’,9,’) are polar co-ordinates relative to O’. 


In general, the interpretation of (57) is somewhat complicated, but the value of 
9%,” is required only on the wing, where 6,’=7. In this case the expressions for 
F’(x) and F’(-*) are particularly simple. It can be seen from (54) that the 
derivative under the integral sign vanishes, leaving 


—ap(x’+r,’) 


F’ (x)= - F’(- x)= -4(x (59) 


Writing € for r,’ in this case, so that € is the distance from O’ along the wing, then 


” rem 60 
1° = o {ap + (5) - (60) 
0 
If the relation 
e~"™du 


K, (p2)= (61) 


is substituted in (60), 9,” becomes 


=- — Nos | x by} 
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—apu 


— hg: | {€+4 (@— b) sin? 6}3 


So 


stew 


The interpretation of this expression is 


z—a(b+€£) 


(s)—h,.~ (s)} ds 


(63) 


When the Mach cone from O’ meets the body, there will, of course be a further 
violation of the boundary conditions there. However, by (63), this will only be 
O(a) and such an error has already been allowed for in the boundary condition 
for ¢,’ (see equation (41)). The solution upstream of the trailing edge is thus 
completed. (The possibility of the cone from O’ running off the opposite tip of the 
wing before reaching the trailing edge, which would produce a significant violation 
of the boundary condition there, is not considered. It is hoped that the present 
analysis will cover most cases of practical interest.) 


There will be a further modification in the anti-symmetrical part of ¢,’, and in 
¢,” and ¢,’, downstream of the trailing edge. The field in this region is determinate 
if it is assumed that the Kutta-Joukowski condition, that the flow leaves the trailing 
edge smoothly, is satisfied. For the present purpose, however, the explicit form is 
not required; for in this region the terms mentioned are not more significant than 
they are upstream of the trailing edge, and it will appear later that these terms do 
not contribute to the forces on the body within the present approximation. 


9. The Axial Force on the Body 


From equation (6), the pressure is given by 


2 Y 


_ 
and the velocity is given by equation (2). On the body ¢,,=R.+O (a’ loga~') and, 
since R. is O(a’) in the region where ¢, and $y are non-zero, it is found that 


2 
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The pressure coefficient on the body is thus 


C,= = . (66) 


and since ¢, is an odd function of 6, it follows that the axial force is given by 


2 
dé {o,.+ S, (z)dz+0O (n?a? + +a* log? a) (67) 


0 


where / is the length of the body (for a body with a flat base there will be a further 
contribution which is not included in this expression). 


(A). Contribution of the terms ¢2+R°/2 
On the body we may write, from the interpretation of (18), 


1 
- Sia (s) log 


2(z-s) 


ds + O (a* log’ a), 


and the contribution to the axial force of these terms is then 


l 
S.(z) a{ S,.(s) log (z s) ds + [s. (z) S..(z) log -tR,’S, dz+ 


0 0 0 
+ log? a) 


t s 
18.0 | S.(2)dz | Su(s)log ds+ 
lfd 2 2 
4S.? log =| dz +O log? a) 
0 
t 


+O(a'log?a). . (69) 
since S,(0)=0 and S,())=O (a’). 


D 
$p,U? | aa RRC,dz 
0 
3) 
er 
be 
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(B). Contribution of the term 4: 


If the inner integral in equation (46) is integrated by parts, it will be seen that, 
even if there are discontinuities in h,., 9,’ is at most O (aloga™*) on the body*. Its 
contribution to the axial force is thus 


0 


*The expression (46) for 9,’ itself will not be reliable near a discontinuity in h,,+. Its 
derivatives are in fact singular along the Mach cone emanating from the discontinuity. The 
same trouble is experienced in the potential for flow past bodies of revolution having a discon- 
tinuity in the meridian cross section. Equation (69), for example, is clearly invalid in this 
case. Lighthill®), however, has shown that, for the body of revolution, a discontinuity 
affects appreciably only a region immediately behind it whose length is O(a). The modified 
result for the drag is that which would be obtained by using the full (linear) relation for 9, 
(equation (24)) instead of the simplified relation (equation (68)) which ceases to be valid. In 
the present problem @,’ does not contribute significantly to the drag even when simplified, 
for the logarithmic singularity occurring is integrable. Lighthill's more detailed analysis, 
however, is an assurance that the contribution to the drag from ¢,’ can safely be ignored within 
the present approximation. 


The contribution of ¢,» is 


and this may be simplified to 


| + + O (na*) (71) 


c 


jp?" 


where the + and - signs refer to the upper and lower surfaces of the wing 
respectively. If, for example, the flow over the wing in the region of the axis is two- 
dimensional, then 


and (71) becomes 
c+1 


ec 


10. The Lift Force on the Body 


The terms 29,. and R.’ in equation (66) for C, are symmetrical about the plane 
y=0, and so the total lift force acting on the body, perpendicular to the axis, will 
be given by 


=2] sin 69 | + Wax) Rdz+O . (74) 
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The term apR?’K, (apr)sin@ (equation (32)) appearing in 9%, is, by (34), O(a) on 
the body. Also, by equations (46) and (48), both ¢,’ and ¢,’ are O(aloga™*). Thus 


Now wy: is zero on y=0 behind the trailing edge, since it is anti-symmetrical 
about this plane, and continuous. Also ¢,~ consists of a symmetrical part, which 
does not contribute to the lift, and an anti-symmetrical part whose z-derivative is 
zero on y=0. Thus the significant term in (75) may be written 
4al (c+ 1)— (C)} (C+ 1) + 


For example, two-dimensional flow past the wing in the vicinity of the axis would 
be described by 


= (z- ay) 


1 
= ay) 
and (76) would then be 
8 
hit C+D} n+ ¥] 


11. The Axial Force on the Wing 


The pressure coefficient on the wing is 
Cy= O (n? + + log? a). 
The force Dy, in the direction of the body axis, is therefore given by 


D 


n.C, drdz 


= -2n | hz Por + 1912 + drdz + O (n° + + nat log?a), . (80) 


the integral being over the upper and lower sides of the plane y=0 occupied by 
that part of the wing outside the body. 
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(A). Contribution of the terms 1%iw+¥ow 


The contribution from 7¢,~+¥%.w. the potential field for flow past the wing 
alone, which is not the concern of the present paper, is excluded from expression 
(80). It is simply noted that their contributions to the drag are O(n’) and O (ny) 
respectively. It should be remembered that the drag of that part of the wing inside 
the body is to be omitted. It is O(y’a+ a) and comparable with terms retained. 


(B). Contribution of the term $v: 
Two particular cases are considered : — 


(i) h, independent of r (implying no sweepback at the trailing edge) and wing 
tip completely outside the Mach cone from the nose of the body. 


The contribution from the upper surface of the wing is then 


e+1 z/a z-ar 


2n S..(s)ds ‘ 


z—aa 


am fn (5. (s)ds| sin~ | +O (ya* log? a), 


c 0 


and since the variation in § is O(a*) downstream of the leading edge, this may be 
written 


c+1 e+1 


| h.*S,(z)dz- 24 | S.,(s)ds h.+ 


c 0 


+ O log?a) 


c+1 c 


h.+S,(z)dz+ - mt S, (s) ds | log (z - s) dh,+ +O (na* log? a) (82) 


c 0 c— 


where the form of the limits of the last integral implies that discontinuities in h.* 
at the leading and trailing edges contribute in its evaluation. The total contribution 
from both upper and lower surfaces is thus given by 


c+1 (e+), 


(2) {het + 74 | log (het + 


c— 
+ O(na* log? a). . . (83) 
(ii) h, independent of r and wing tip completely inside the Mach cone from the 
nose of the body. Let the wing be rectangular and of semi-span b measured from 
the axis. Then the contribution from the upper surface of the wing is 


zla z—ar 
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and, by ectiied with (81) and (82), this must be 
c+1 C+D4 


(2)de+ 2% 14 (s)ds | log - s)dh,* - het, (z ba) dz+ 
c 
ce+1 


2 b 
| h.*dz S..(s) sin~" (=) ds+O(na‘log?a). . (84) 


ec 
The total contribution is thus given by 


c+1 
(z) - S,(z- ba)} +h,-} dz+ 


(s)ds { log (z-s)d {h.* + 


0 c= 


c+1 z—ba 


+ (h.+ dz S..(s) sin=" ds + O (na* log? a). 


(C). Contribution of the terms and 


(i) h, independent of r and wing tip completely outside the Mach cone from the 
junction of the wing leading edge and the axis of the body. 


The contribution from the upper surface of the wing is 


z-—c 
c+1 @ 


41? | h.tdz 
c a 


Now, from the interpretation of equation (50), 


z-—ar 


_ 2a h,* (s) ds 


except in a region near the body whose width is O(a), where 9,’ is O(aloga™'). 
This will not contribute to the drag to within O (7’a? log a~') so that (86) becomes 


dr - + O (n?a* log 
0 


c= 
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we (s) i 


C= 0 


Pry . (88) 


(h.+)? dz +O (n?a? loga-'). ‘ . (89) 


The total contribution from the upper and lower surfaces is therefore given by 
c+1 


(ii) h, independent of r and the semi-span < 2~' (Fig. 2). The wing is assumed 
to be rectangular and of semi-span b downstream of the point at which the Mach 
cone from the junction of the wing leading edge and the axis runs off the tip of 
the wing. 


The contribution of the upper surface is now 


z—c z-—c z-c—b 
c+1 c+1 c+1 


| +dz | hetdz | 9, ~ | | +0 loga-) 
b 


c e+ab 


e+ab 
(91) 
where 9,’ is given by (87) and ¢,” by (63). 


The first term has already been evaluated. The second, by analogy with (88). 
can be written 


e+1 


Bay" j hitdz j dh.+ (s) 


{(z- 


b 
e+1 z—ab 
“er. h,* (z)h.* (z- ab) dz - 8 | h,+dz | sin-! (s). 
(92) 
Finally the third term is equivalent to 


z-—c 
c+1 


d {h,* 
h dé J {(z- - a? (b— 


(93) 


c+ab 


308 


c+1 
a 
c 
c 
; T 
( 
a 
a 
Pp 
r 
t 
I 
|_| 


SUPERSONIC WING-BODY COMBINATIONS 


Collecting the terms for the upper and lower surfaces, and remembering that 9,” is 
anti-symmetrical about the plane y=0, gives finally the expression 


c+1 


(z) {h.+ (z)—h.* (z - ab)} (2) (2)— (z— 2b)} 


e+1 z—ab 


n? b 
| (h.+ +h.-) dz | sin-' () - 


c+ab 


4an? 
a thet (9) + 
e+ab c= 


+O (n?a* loga™'). . (94) 
There is now only the contribution from ¢, to be considered. The term 
xpR*K , (apr) sin @ occurring in ¢, is, of course, zero on the wing. Furthermore, by 
(50), 9.’ is O(a? log a~') on the wing except for a region near the axis whose width 
is O(a), where 9,’ is O(aloga~'). The contribution to the drag from ¢.’ is thus 
at most O loga~'). 


12. The Lift Force on the Wing 


As for the drag, the contributions from ¢,, and ¢,, should be restricted to that 
part of the wing outside the body. These contributions will be O(n) and O(y) 
respectively. 


Since 9, is symmetrical about y=0, and since the contribution to the lift from 
,, is at most O (a* log a~'), the only further contribution to be considered is that 
from 9,’+9,”. This will be 


2n +912") drdz 


taken over the two surfaces of the wing. The error will be O(n? +¥? +a‘ log’ a). 


(i) If the wing tip is outside the relevant Mach cone and the trailing edge is 
perpendicular to the axis of the body, the contribution from the upper surface is 


z—c 


c+1 
This compares with equation (86) apart from a change of sign and the suppression 
of nh,*. Thus, by analogy with (89), this must be equivalent to 
c+1 
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Since ht (z)= — h~ (z) at z=c,c+1, the total lift contribution is given by 


8 2 -1 


(ii) If the wing tip is partly inside the Mach cone and parallel to the main 
stream there, the contribution from the upper surface becomes 


z-c—b 
c+1 c+1 


—4y (a nie dr + 4 | d z | . (97) 


e+ab b 0 


where the first term is given by equation (95), and the second, by analogy with (92), 
is equivalent to 


z—ab c+1 z—ab 
_ Bay + = ( 2b 
dc h,* (s) cos ds 
c 


c+1—ab 


c 


z=8+4ab 


e+1—ab 


8an 
h,* (z)cos-! 


c 


Finally the third term may be compared with (93) and is equivalent to 


e+1 


ae 2 {2 1} Ja ther 


e+ab c— 


e+1—ab 


c 


Collecting the contributions from the upper and lower surfaces gives the following 
expression for the lift 


;= [ 21° (c+ 1)—ht (c+1-2b)+h- (c+1—2b) | + 


h,~) sin 1 ¢ dz+O(na* loga™'). 
¢ 
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13. Comments 


The calculation of the error terms at each stage, which a comparison of the 
orders of the terms retained and the error terms shows to be highly desirable, has 
resulted in a somewhat involved analysis, although the practical application of the 
final results should not be difficult. 


The only contribution to the lift from the body, for example, arises from the 
potential of the wing alone (equation (76)) and for two-dimensional flow past the 
wing in the neighbourhood of the axis this reduces to the simple relation (78), 
namely 


where the first term inside the square brackets is simply the wing-body setting. 


This is just the contribution which would be expected from the part of the 
wing lying inside the body for flow past the wing alone“. If this is a reliable guide 
in all practical applications then this term may be absorbed into the wing lift, if 
this lift is calculated on the basis of the complete wing including the portion lying 
inside the body. If the wing is symmetrical about the plane y=0, there is no 
further contribution to the total lift within the approximation. For an asymmetrical 
wing, or a symmetrical wing at incidence relative to the body axis, there is, in 
addition, the contribution from the interference potential on the wing. If the aspect 
ratio is greater than 2(M?— 1), then this is given by equation (96), namely 


4p,U*an 
(c-+1)-ht 


It will be seen that this is just sufficient to annul the contribution from the body lift 
arising from asymmetry (see equation (78) ). 


Speaking generally, therefore, if both the wing and body are at the same 
incidence relative to the main stream, the lift is effectively that of the whole wing 
when isolated from the body. If the incidence of the wing is further increased 
relative to the axis of the body, the increase in lift is that associated with that part 
of the wing outside the body, for flow past the wing alone. It is noteworthy that 
this loss of lift does not arise directly from the shielding effect of the body on the 
wing, for there is a contribution to the lift from the body and this contribution is 
equivalent to the lift of that part of the wing inside it. The decrement arises from 
a loss of lift on that part of the wing outside the body and inside the Mach cone 
from the junction of the wing leading edge and the body axis. 


The error terms are all included in O(n? +¥* + a@¢ log? a). 


The drag force is more interestinz. Apart from the contribution of the wing 
alone, there is the contribution of the body due to its own potential field, and given 
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by equation (69). The result is well-known®:'’:'” and requires no comment save 
possibly on the occurrence of the first term in (69). This is not new and can be 
deduced from Ward’s expression for the drag of slender bodies“”. This 
contribution is O(a’). 


In addition there is the further contribution arising from the effect of the 
potential field of the wing on the body given by (71). It may well be that the 
simplified expression given in equation (73) is sufficiently accurate for most practical 
purposes, for this term is O(ya*) and in any case S, is likely to be small in this 
region. 


From the wing there is a contribution due to the potential field for the flow 
past the body. If the wing tip lies outside the Mach cone from the nose and the 
wing ordinate is independent of r, this is given by equation (83) and is O(ya’*). If 
the wing tip lies inside the Mach cone and the wing is rectangular, then equation 
(85) is used; this theoretically is O(na?) and could thus compare with the 
contribution from the wing potential itself if, for example, y,=O(a’). It may seem 
surprising at first that the effect of the nose is less significant when the aspect ratio 
is large, than when it is small. The reason is that the integrated pressure over an 
elementary spanwise strip of the wing, arising from the perturbation field of the 
nose upstream of the leading edge, is independent of z provided the wing ordinate 
is independent of r and the wing tip lies outside the Mach cone from the nose of 
the body. There is thus no contribution to the drag, since a constant force acting 
round a closed contour, and normal to it, has a zero resultant. This result, however, 
does not hold if the strip encounters the wing tip before the outer boundary of the 
Mach cone (or cuts through the junction of the wing and body—hence the second 
term in equation (82)). 


If §.(z) is non-zero only for a short distance near z=0, then the contribution 
to the drag from Dy,, given by (85), can be made small by fixing the wing well 
downstream of the nose, as would be expected. For then the expression 
{S.(z)-S.(z—b)} in the integrand of the first term is zero throughout the range of 
integration. In the second term we have 


If the fact that h+ =—h~ at the leading and trailing edges is used, it will be seen 
that the integral is zero and the original expression is thus O(1/c’). 


In the third integral of (85), sin- {zb/(z-s)} is effectively independent of s 
where S,,(s) is not zero. The integral of S,,(s) then vanishes. 


Finally, there is the contribution D,, from the interference potentials 9,’ and 
¢,”, which is O(n?a). If Dy, is small and the body is cylindrical downstream of 
the leading edge, then this is the only significant extra term, apart from the drags of 
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the wing and body due to their own potential fields. An interesting point is that, 
at least for wings of aspect ratio greater than 2(M*-1)-}, it is always negative 
(see equation (90)). Thus a decrease in drag may be expected, quite apart from 
the shielding effect of the body on that part of the wing inside it. It should be 
remembered, however, that if the wing is well forward near the nose, the effect of 
D., May turn the scales completely. 


The similarity of equation (90) with Ackeret’s result”? for the drag of a two- 
dimensional wing is striking. In fact, if the body consists of an infinite cylinder of 
radius a, and the wing is a portion of a two-dimensional wing for which Ackeret’s 
result can be used, then the drag of the complete wing alone will be 


c+1 


a 


Dwing = * 


where 25 is the span. 
The drag of the wing-body combination, however, is by this result and by (90) 


Dev=  (b-2a) | dz 


Daw 


The decrement arises from the shielding effect of the body on that part of the 
wing inside it, and from the interference between wing and body. It is equally 
divided between them. 


The error terms in the drag are all included in the expression 


if b>(M?-1)-}. 


O (n° + nya’ log a~' + log’ a + ny? + Va? +a° log’ a). 
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The EDITORIAL BOARD has much pleasure in announcing 
that, beginning with Volume V, The Aeronautical Quarterly is again 
to be published four times a year. 


Upon the interval between issues largely depends the time lag 
between submission of a paper and its publication, and the Board 
hopes that by reverting to quarterly publication, this time lag will be 
reduced correspondingly. The interval between submission and 
acceptance (or rejection) of a paper should also be much reduced. 


Many of the papers submitted in the past have been mathematical 
in character and the Board considers that there is a tendency to include 
in the text of papers over-much mathematics; to facilitate the reading 
of papers they would prefer, where possible, to see mathematical 
derivations of any length included in an Appendix. 


In certain fields the number of papers submitted has been satis- 
factory but in others, few papers have been received. The Editorial 
Board feel that there are many sources from which original papers 
on a wide variety of subjects should be available and they wish to 
emphasise that the Editors will welcome original papers on all 
branches of the aeronautical and allied sciences. 


The Board wish to encourage especially, more papers describing 
iii 
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experimental work and techniques, both to give The Quarterly a better 
balance and to attract a wider circle of readers. Experimental workers 
in Universities, in Research Establishments and in the Industry could 
find The Quarterly a valuable publication. 


The Editorial Board 


